domain: phys sde

text:
Faustine Mascaut
email: faustine.mascaut@univ-lille.fr

Olivier Pujol
email: olivier.pujol@univ-lille.fr

Bert Verreyken
email: bert.verreyken@aeronomie.be

Raphaël Peroni
email: raphael.peroni@univ-lille.fr

Jean Marc Metzger
email: jean-marc.metzger@univ-reunion.fr

Luc Blarel
email: luc.blarel@univ-lille.fr

Thierry Podvin
email: thierry.podvin@univ-lille.fr

Philippe Goloub
email: philippe.goloub@univ-lille.fr

Karine Sellegri
email: k.sellegri@opgc.univ-bpclermont.fr

Troy Thornberry
email: troy.thornberry@noaa.gov

Valentin Duflot

Pierre Tulet

Jérôme Brioude

Aerosol characterization in an oceanic context around Reunion Island (AEROMARINE field campaign)

Keywords:

Optical properties, vertical distribution and transport pathways, from the marine boundary layer to the free troposphere, of marine aerosols in a pristine environment are examined.

• Aerosol size does not exceed the accumulation mode.

• A sketch is proposed as a characterization of marine aerosols distribution. Oceanic and insular influences in the aerosol content are separated.

• It is argued that the AERONET station at St Denis (Reunion Island) is well representative of marine conditions.

• With the thermodynamics given by a microwave radiometer, the results will be useful for water vapor-aerosol-cloud physical processes modeling in a pristine ocean.

Introduction 19

Because of their direct and indirect radiative forcings, 20 atmospheric aerosols have a major impact on the climate.

21

These forcings are still poorly understood and lead to uncer-22 valentin.duflot@univ-reunion.fr (V. Duflot); pierre.tulet@aero.obs-mip.fr (P. Tulet); jerome.brioude@univ-reunion.fr (J. Brioude) ORCID(s):

tainties that have persisted in models since the 1990s [START_REF] Myhre | Anthropogenic and Natural Radiative Forcing[END_REF]. One of the largest uncertainties of the aerosolcloud system is the background concentration of natural aerosols, especially over clean marine regions [START_REF] Andreae | Aerosol-cloud-precipitation interac-847 tions. part 1. the nature and sources of cloud-active aerosols[END_REF]. Oceans cover about 70% of the Earth's surface and are an important reservoir of marine aerosols (mainly sea salt and organic aerosols). Sea Salt Aerosols (SSA) are one of the largest contributors to global aerosol loading and therefore they play an important role in global climate. Also, they were proposed to be a major component of primary marine aerosol mass over the regions where wind speeds are high and/or other aerosol sources are weak (Gantt and Meskhidze, 2013;[START_REF] O'dowd | Marine aerosol production: a review of 993 the current knowledge[END_REF]. [START_REF] Luo | On factors controlling marine boundary layer aerosol optical depth[END_REF] showed that the at 532 nm of SSA in the Marine Boundary Layer (MBL) is governed by different physical factors: the surface wind speed, Sea Surface Temperature (SST), MBL height, lower troposphere stability, and relative humidity (impact on size and optical properties of sea salt particles). However, there is a notable lack of data, in particular over the oceans and in the southern hemisphere, describing the characteristics of aerosols such as optical properties, size distribution, temporal and spatial variabilities... [START_REF] Ramachandran | Spectral aerosol optical characteristics during the 1006 northeast monsoon over the arabian sea and the tropical indian ocean: 1007 1. aerosol optical depths and their variabilities[END_REF]. [START_REF] Pant | The concentration and number 997 size distribution measurements of the marine boundary layer aerosols 998 over the indian ocean[END_REF] carried out measurements of the total number concentration and the size distribution of aerosols over the Indian Ocean in 2004. They observe that the aerosol concentration-wind speed correlation coefficient depends on the latitude and has a maximum value where the winds are the strongest.

The southwestern Indian Ocean has been identified as a region with frequent pristine conditions (where land and 52 human activities have little impact) that can reasonably be 53 considered to be close to the preindustrial conditions. Few 54 data have been collected in this region [START_REF] Pant | The concentration and number 997 size distribution measurements of the marine boundary layer aerosols 998 over the indian ocean[END_REF].

55

However, it is a crucial reference point to quantify the back-56 ground concentration of natural aerosols and the contribu- 68 Indeed, being in an oceanic environment and far from con-69 tinents, the island is in a strategic location for carrying out 70 measurements in a clean region, and also for the validation 71 of spatial measurements. In addition, the Maïdo Observa-72 tory (located at 2.2 km above sea level (a.s.l)1 ) allows: (i) to 73 take measurements directly in the free troposphere at night 74 (Guilpart et al., 2017;Foucart et al., 2018) and (ii) to perform In this text, wavelengths are given with respect to the vacuum dependence (Karol et al., 2013). The Angström exponent is determined from the law of Angström:

1 ∕AOD 2 =
( 1 ∕ 2 ) -where is the aerosol optical depth at the wavelength (AngstrÖm, 1961). This exponent describes how the varies with and so provides information on the size distribution of the aerosols (Kusmierczyk-Michulec et al., 2002). sphere and to investigate fruitfully a wide variety of weather phenomena related to water vapour [START_REF] Louf | Seasonal and diurnal water vapour distribution in the sahelian area from microwave radiometric profiling observations[END_REF]. A zenith-looking infrared ceilometer provides, together with the temperature profile retrieved from the MWRP, an estimate of the cloud-base height. The MWRP is also equipped with in situ sensors for ground level measurement of temperature, water vapour and pressure [START_REF] Louf | Seasonal and diurnal water vapour distribution in the sahelian area from microwave radiometric profiling observations[END_REF]. It was on the university campus between December 12th, 2018 and March 11th, 2019.

Particle counters

Models and reanalyses

AROME

AROME-Indian Ocean (Bousquet et al., 2020) is used in this study in order to obtain the horizontal wind fields at different altitudes, at the places and dates where the flights were performed. This model is an adaptation of Météo-France's operational model AROME [START_REF] Seity | The arome-france convective-scale operational 1021 model[END_REF] PART has been coupled to the Eulerian models AROME (Verreyken et al., 2019). In this paper, we also use a version of FLEXPART that has been coupled to meteorological output from Meso-NH.

Meso-NH

CAMS reanalyses

Copernicus Atmosphere Monitoring Service (CAMS, previously MACC, https://atmosphere.copernicus.eu/) is an atmospheric model that simulates the mixing ratios of various aerosols, or thermodynamic parameters (for example humidity, wind, temperature ...) on a large scale and regional scale. The data assimilation system used for CAMS is based In our study, we use the mixing ratios of these different species [START_REF] Mallet | Marine aerosol distribution and variability over the pristine southern indian ocean[END_REF] as well as the . Mixing ratios and aerosol concentrations are directly related, allowing comparison between in-situ measurements and CAMS data.

Results

From global to local scale

≈ 8.314 J K -1 mol -1 . Since = ∕ , it ensues directly ∕ = 3 ∕(4 3
).

The results of this conversion are summarized in Table 1. These orders of magnitude are realistic values. Another We note for these flights that the SSA are not located on 321 the same side of Reunion Island depending on the day. We 322 will see later that this is explained by the wind regimes. Table 1 Range of concentrations measured by POPS, TSI and MAGIC and range of SSA concentrations (size ranges: 0.03-0.5 m; 0.5 -5

m and 5 -20 m) from CAMS reanalyses in the pixel where the flight took place, for the six flight dates.

Date MAGIC (#/cm 3 ) TSI (#/cm 3 ) POPS (#/cm 3 ) SSA 0.03-0.5 (#/cm 3 ) SSA 0.5-5 (#/cm 3 ) SSA 5-20 (#/cm 3 ) 03/13 200 to 2.10 3 100 to 10 3 5 to 100 0.3 to 1.4.10 3 1.2.10 -3 to 0.1 4.6.

337

From the radiometer measurements, we observe that rel- indicate a depletion of small particles in favor to larger ones.

Aerosol concentration measurements

Figure 6 presents the aerosol concentration measurements by MAGIC (particle size larger than 2.5 nm), TSI (particle size larger than 10 nm) and POPS ( particle size larger than 132 nm). POPS concentration profiles show that the sampled particles (accumulation and coarse modes) have a maximum number concentrations of around 10 2 cm -3 between 0.5 and This means that small particles (nucleation mode) are twice (in concentration) than larger particles (Aïtken mode) below 1.5 km while above, there is no small particles (MAGIC-TSI ratio close to one). Such values of number concentration in the nucleation mode, between 10 3 and 10 4 cm -3 , have been measured also in the Mediterranean (Eleftheriadis et al., 2006).

For the flights in group I, the concentration vertical profiles obtained from MAGIC and TSI present similar shapes.

On average, MAGIC concentrations are 2 to 3 times higher than TSI concentrations. The highest concentrations are found in the marine boundary layer below 0.7 km a.s.l for flight F1, and below 0.5 km a.s.l for flight F2 and F3, with values ranging from 2000 to 10000 cm -3 . This difference of 200 m in the marine boundary layer may be due to the fact that flight F1 occurred one hour later in the morning compared to flights F2 and F3, so the boundary layer may have time to develop a bit more.

473

For flight F1, the MAGIC and TSI concentration profiles 474 are rather constant at 2.10 2 cm -3 above 1 km. For flight F2, 475 the concentrations are constant above 1.5 km at 4.10 2 cm -3 .

476

For flight F3, a value of 4.10 2 cm -3 is found above 2 km, 477 and a value of 10 3 cm -3 is found in the layer between 0.8 ing on the flight, these modes have a lower concentration and larger particles (size greater than 132 nm -accumulation mode) are dominant. The exception in the marine boundary layer for flight F5, i.e. below 0.5 km, means that this layer is poorer in terms of small aerosols.

Origin of the air masses

To further evaluate the origin of the air masses identified in the previous section, we used FLEXPART mesoscale backtrajectories to distinguish air masses influenced by local terrestrial emissions from those with marine origin or representative of the regional background. Here, we have chosen a backtime of 12 h because it corresponds to the domain of AROME.

The tables A1 and A2 in the appendix presents statistical products from the FLEXPART output.

The scientific meaning of the FLEXPART outputs needs to be briefly reminded in order to avoid misunderstandings about the result presented in those tables. When we release an air mass at a given altitude , FLEXPART is able to trace back its probable trajectories, over a geographical grid, during an user-decided time interval. From that, it is then possible to identify (and reckon) the oceanic pixels, viz. the trajectory grid-points located over the ocean, and similarly for the island pixels which are trajectory grid-points located over Reunion Island. In other words, we can quantify how much could the ocean contribute to, or impact, the aerosol content of the considered air mass.

Beside the geographic origin, it is possible to determine the tropospheric layer from where an air mass originates.

530

The results from the previous sections helped us identify 531 three layers: (L1) below 0.5 km, (L2) between 0.5 and 1.5 km, 532 and (L3) above 1.5 km. For instance, in Table A1, the results

533

after 12 h of backward simulations for an altitude release at A1. It is also at this altitude that PLASMA measures the most important (during the descending phase).

We are summarized with typical orders of magnitude the results given in Tables A1 andA2. These conclusions are summarized on Figure 7. The vertical layering we found, and in particular the predominance of the oceanic aerosols, are explained by the interaction between the wind field and The difference in concentration between MAGIC and POPS is much larger over St Denis than over the ocean. This means that, as indicated by the PLASMA optical measurements, small particles (size lower than 132 nm) are the dominant mode over St Denis. The MAGIC -TSI differences in concentration further indicates that the smallest particles are those that dominate. In contrary, the MAGIC -TSI difference is close to zero and confirm that larger particles are found above 1.5 km a.s.l over the ocean north east of the island.

The wind direction given by AROME (not shown here) indicate (i) at 1 km (overflight of St Denis) a southeast wind of around 10 m s -1 coming from the island and (ii) at 1.5 km, when the plane begins these spirals above the ocean, a southeast wind from the ocean with a speed of about 8 m s -1 .

The FLEXPART backtrajectory results (table A2 in the appendix section) indicate that for a release at 1.5 km, the air The time series (POPS, TSI and MAGIC measurements) dur-ing flight in the free troposphere are of the same order of magnitude as the Scanotron night measurements made at the Maïdo Observatory, in the free troposphere (between 1×10 2 and 7 × 10 2 cm -3 ). The two TSI CPCs (at the Maïdo Observatory and on the plane) also measure identical concentrations above 2 km (around 10 2 cm -3 ).

We can conclude that the nighttime measurements at the Maïdo Observatory (at 2.2 km a.s.l) are representative of the daytime measurements (during flights) in the free troposphere and allow sampling of purely marine aerosols.

Conclusion

In

57

  tion of natural emissions to the changing climate. In pristine 58 regions, SSA are dominant and concentrations are relatively 59 low (e.g.[START_REF] Mallet | Marine aerosol distribution and variability over the pristine southern indian ocean[END_REF]. 60 Reunion Island in the southwestern Indian Ocean can be 61 considered as a background aerosol pristine environment un-62 der trade wind conditions (Koren et al., 2014), mostly dur-63 ing the wet season from December to April. For details 64 about wind circulation in the Southern Indian Ocean, includ-65 ing Reunion Island, see the statistical study by Mallet et al. 66 (2018) and refs therein. Reunion Island is also a unique site 67 in the southern hemisphere for making aerosol observations.

  Photomètre Léger Aéroporté pour la Surveillance des Masses

  on the ECMWF' Integrated Forecast System (IFS). Satellite observations are implemented in this model and allow the study of the atmospheric composition (chemically reactive gases, aerosols, greenhouse gases) at global scale (Morcrette et al., 2009). CAMS allows in particular to differentiate various types of aerosols such as: (i) Sea Salt Aerosols (0.03-0.5 m; 0.5-5 m and 5-20 m) and Dust aerosols (0.03-0.55 m; 0.55 -5 m and 5 -20 m) divided into three size ranges, (ii) Black Carbon and Organic Matter divided into two modes (hydrophobic and hydrophylic) and (iii) Sulfate.

291Figure 2

 2 Figure2shows the CAMS mixing ratios (for three

Figure 1 :

 1 Figure 1: CAMS reanalysis (0.25°/0.25°): mixing ratio of SSA 0.03-20 m (top left), sulfate (top right), Organic Matter and Black Carbon (bottom left) and Dust Aerosols 0.03-20 m (bottom right) at 850hPa. Date: 03/22/2019 at 06∶00∶00 UTC.

  precise approch would be to consider a size distribution if it were fully available. The POPS, TSI and MAGIC size ranges of measurements (respectively 132 nm-3 , < 10 nm and < 2 nm) correspond to the two smallest size ranges of the CAMS reanalyses (SSA between 0.03 m and 5 m). The CAMS concentrations retrieved (SSA 0.03-0.5 and SSA 0.5-5 ) are of the same order of magnitude as the concentrations measured by POPS, TSI and MAGIC (see Tab.1). The largest SSA concentrations (SSA 5-20 ) are negligible compared to 319 the other two CAMS size ranges.

  320

Figure 3 Figure 2 :

 32 Figure 3 present the averaged profiles (for March 2019)

Figure 3 :Figure 5

 35 Figure 3: Average profiles of temperature and relative humidity from microwave radiometer and from CAMS reanalysis (20.8°S; 55.2°E) for March 2019.

Figure 4 :

 4 Figure 4: Flight trajectories and altitudes on 03/13/2019, 03/15/2019, 03/27/2019 (Group I) and 04/11/2019, 04/18/2019 (Group II). The sixth flight is a particular case (see text for details).

Figure 5 :

 5 Figure 5: and Angström exponent ( ) vertical profiles (PLASMA measurements) for flights F1, F2, F3 (Group I) and flights F4, F5 (Group II).

  , and around 10 to 3 × 10 2 cm -3 above. It confirms 425 that larger particles are found between 0.5 and 1.5 km in al-426 titude, in agreement with the conclusions deduced from the 427 optical measurements. Below 0.5 km a.s.l, the POPS con-428 centrations increase with altitude during the ascent phase.429The shape of the measured concentration profile can show 430 a maximum (group I) or a minimum (group II) in the layer

Figure 6 :

 6 Figure 6: Particle concentration profiles (from POPS,TSI3007 and MAGIC200 instruments) during F1, F2, F3 (Group I) and F4, F5 (Group II).

  , only MAGIC was available. Similar fea-480 tures can be found in the profiles of flights F4 and F5 com-481 pared to group I. Furthermore, contrary to the aerosol con-482 centration measured by POPS, the aerosol concentrations are 483 rather similar between the ascent and descent phases. 484 Other important information are obtained when compar-485 ing the POPS and MAGIC profiles. Indeed, the concen-486 tration profile differences "MAGIC -POPS" (not shown) 487 present relatively smaller values above 0.5 km. More sig-488 nificant differences are obtained below this altitude. An ex-489 ception occurs for flight F1 since the critical altitude is 1 km 490 because of a more developed boundary layer, as already sug-491 gested above. Another exception is for flight F5 where sig-492 nificant differences are between 0.5 and 1.5 km and small 493 differences elsewhere. For group II, no significant differ-494 ences are noted between the ascent and the descent i.e. same 495 shape and order of magnitudes in the concentration profile 496 differences. 497 In summary, the measurements indicate that the bound-498 ary layer below 0.5 km -0.7 km is much richer in aerosols 499 in the nucleation or Aïtken modes (size lower than 132 nm). 500 Above the boundary layer, and up to 1.5 km to 2 km depend-Aerosol characterization in an oceanic context around Reunion Island (AEROMARINE field campaign)

534

  200m for flight F1 show that 88% (first row, last column) of 535 the air mass originated from grid cells over the ocean, and 536 12% from grid cells over Reunion island. In addition, the ori-537 gin of the air mass can be analysed in terms of geographic 538 location and height. Hence, 43% of the air mass originated 539 from 0-500m in altitude (layer L1), 38% from 500 to 1500m 540 in altitude (layer L2), and 7% above 1500m (layer L3) and 541 above the ocean. Therefore, the origin of this air mass is 542 mainly oceanic. 543 We have thus used FLEXPART to calculate the origin 544 of air masses from different altitudes. For group I, we can 545 see that the origin of the air masses is mainly oceanic since 546 more than 73% of the backtrajectories originate from grid 547 cells above the ocean. 548 In contrast, 28% of the air mass was located above Re-549 union island, and within the layer L2 and L3. For group II, 550 the results are different since the air masses present a sig-551 nificant origin above Reunion island. In particular, the air 552 masses released at the 200 m and 1 km altitudes for flight F4 553 have a dominant origin above Reunion island. 554 Comparing the altitude where the backtrajectories are re-555 leased, and the distribution in the vertical of the backtrajec-556 tories give additional information on the vertical transport 557between layers. Let us first look at oceanic pixels for group I. First, we note that the air masses with a dominant marine origin stayed, for a large part, within the same layer they originated from. For example, for a release at 200 m, 43% of the air mass with marine origin were in the layer L1. Excepted for some rare exceptions (underlined in tables), the layer that presents the highest percentages (in bold fonts in tables) includes the altitude of release. However, exchanges of air masses (in italic fonts in tables) between layers are not negligible, even if they are not dominant. For instance (flight F1), for a release at 1 km, 29% of the air mass with marine origin come from the layer L1. This means that the aerosol content in the lower troposphere is also impacted by the mixing between contiguous layers. This result holds for the air masses originating from Reunion island, although the terrestrial component in the lower troposphere is mostly significant in the layer L2 (percentages greater than 10%). For the group II flights, we have the same results although we recognize that the insular influence is more present, especially at low level for flight F4. This may be correlated with MAGIC and PLASMA measurements. Indeed, for flight F4, the maximum MAGIC concentrations are measured below 0.5 km a.s.l, which corresponds to a majority of island pixels in the Table

588Figure 7 :

 7 Figure 7: Descriptive diagram of aerosols inputs and the exchanges between the atmospheric layers L1 (below 0.5 km), L2 (0.5 -1.5 km) and L3 (above 1.5 km).

Figure 8 :

 8 Figure 8: Particle concentration profiles from POPS, TSI3007 and MAGIC200 instruments (left), and Angström exponent ( 440∕870 ) profile from PLASMA sun-photometer (right) during flight on 03/22/2019 (F6).

621

  photometer (St Denis) for each flight date

Figure 9 :

 9 Figure 9: Mean , mean and number of data points from AERONET station in St Denis (Level 1.5) for the period from March 17 to May 17 (2010 -2020).

  Figure 10 presents the daily averages of the CAMS

Figure 10 :

 10 Figure 10: measurements from AERONET sun-photometer (Level 1.5) and from CAMS reanalysis (20.8°S; 55.2°E) for March 2019 and April 2019.

690(

  Figure 11 display the mensual averages for March and

  699

761 ••

 761 this paper, we have presented the AEROMARINE field campaign which took place between February and April 2019 off the coast of Reunion Island (southwestern Indian Ocean). This area, identified as a pristine region, is of major interest for the study of marine aerosols, their vertical distribution and their optical properties. During this campaign, a MWRP was deployed in St Denis (93 m in the north of the Island) between mid-December 2018 and mid-March 2019. This made it possible to determine that, during the austral summer in this region, the thermodynamic situation (humidity, temperature, and height of the boundary layer) is relatively stable. In addition, six instrumented flights allowed the aerosols to be sampled from an altitude of 100 m up to 4 km by spiraling above the ocean thanks to an instrumental synergy. The optical properties of the aerosols were measured by the PLASMA photometer and three particle counters (POPS, TSI and MAGIC) measured the aerosol concentrations for 748 different size ranges (accumulation, coarse and Aïtken modes). 749 POPS analysis indicates that almost all of the particles are 750 in the accumulation mode, centered around a particle size of 751 indicate that the aerosols are in the accumu-756 lation and coarse modes, and mainly below 2 km of altitude.757 The FLEXPART simulations enabled to determine the 758 most probable origin of the aerosols measured during the 759 flights. As a result, the aerosols follow the following vertical 760 distribution: Above an altitude of 1.5 km, the sampled aerosols are 762 not substantially impacted by the surface (layer L1 has 763 a relatively little contribution). This is interesting since 764 it allows to quantify the background aerosol concen-765 tration. For all the flights, we have estimated that the 766 number concentrations (in cm -3 ) are 300 (MAGIC), Below 0.5 km (in the MBL), aerosols come essentially 770 from the surface. The origin can be oceanic (33%) or 771 insular (8%). Insular influence are nonetheless due to 772 special events depending on the wind regime (e.g. , 773 the Cap La Houssaye may bring sometimes dust aerosols).774 • The intermediate layer, i.e. between 1.5 and 0.5 km a.s.l, is a layer of mixture: aerosols are mixed with those coming from the lower or upper atmospheric layers. These results meet the initial objectives of the AERO-MARINE campaign: (i) to characterize marine aerosol optical properties and their vertical distribution and (ii) to examine the transport pathways of marine aerosols from the MBL to the free troposphere. It is worth mentioning that the flights were carried out between 04∶00 and 07∶00 UTC, viz. during the transient convection regime between nocturnal and diurnal conditions. Further observational studies and field campaign may be necessary to examine aerosol distributions during purely diurnal and nocturnal regimes, i.e. for well established regimes. The AEROMARINE campaign presented here is interesting in the sense that it documents a transient regime, namely a more complex regime in terms of thermodynamics compared to established ones. The measurements taken during the flights were compared with the CAMS reanalyses. They showed that, like inflight measurements, SSAs are predominant around Reunion Island and that aerosols are mainly located below 2 km. It was also shown that CAMS overestimates the (from 0.01 to 0.09) in this region in agreement with results from Mallet et al. (2018). In addition, a comparison between PLASMA measurements (on the runway) and the AERONET sun-photometer (located in St Denis) as well as a study on the impact of 2020lockdown due to the Covid-19 pandemic on AERONET measurements were carried out.

photometer we used is located on 181 the university campus, at St Denis located in the north of the 182 island. 183 2.2.2. Microwave radiometric profiler (MWRP) 184

	167
	168	2.2. Ground-based instruments
	169	2.2.1. AERONET stations
	170	The AERONET collaboration provides globally distributed
	171	observations of spectral	, inversion products, and pre-
		cipitable water in diverse aerosol regimes. Aerosol opti-
	177
	178	For comparison with PLASMA measurements, only Level
	179	2.0 data quality for	and the Angström exponent ( 440∕870 )
	180	are used. The 2020 data from the AERONET station (St De-
	185	nis) is Level 1.5. The sun-The microwave profiler RPG-HATPRO G5 gives us mea-
	186	surements of the microwave radiation emitted by the tro-
	187	posphere which provides tropospheric vertical profiles (0 -
	188	10 km) of absolute humidity and temperature, with a special
	189	focus on the MBL. It allows to monitor with a high tempo-
	190	ral resolution (1 min) the thermodynamic state of the atmo-

The Portable Optical Particle Counter (POPS) is a 900 g in situ instrument designed by US laboratories NOAA and CIRES that provides aerosol number size distribution (in the size range 132 nm -3 m) using single-particle light scattering

(Gao et al., 2016)

. POPS is a prototype made by a 3D printer to reduce weight. It flew on board a light plane during the AEROMARINE intensive field campaign within the MBL and the free troposphere.

Two Condensable Particle Counters (CPCs) (accuracy:

± 20 %) are used simultaneously to measure the total concentration of particles larger than 2 nm (CPC-MAGIC200) and particles larger than 10 nm (CPC TSI model 3007). The difference of concentration between the two CPCs gives the particle concentration in the size range 2 -10 nm, which is indicative of the recent formation of nanometric particle, i.e. nucleation. The combination of two CPCs to investigate nucleation was proven to be adequate in past airborne studies

(Crumeyrolle et al., 2010)

. Additionally, in synergy with the POPS, the CPC TSI 3007 concentration enables to measure the aerosol concentration in the 10-150 nm size range, which is indicative of grown nucleated particles, and fine 164 marine primary aerosols (size range 50 -100 nm) that dom-165 inate the primary marine aerosol size distribution

(Schwier 166 et al., 2017)

. 172 cal properties are measured at multiple wavelengths ranging 173 from the UV to shortwave infrared. data (accuracy: 174 ± 0.02) are computed for three data quality levels: Level 175 1.0 (un-screened), Level 1.5 (cloud-screened), and Level 2.0 176 (cloud screened and quality-assured) (Holben et al., 1998).

Table 4

 4

	Differences in daily	averages of
	the CAMS model (20.8°S; 55.2°E) and
	of the AERONET sun-photometer (St
	Denis) for each flight
	Flight	-
	F1	+ 0.03
	F2	-0.01
	F3	+ 0.04
	F4	+ 0.09
	F5	+ 0.01
	F6	+ 0.06
	(2018) who determined that the	given by CAMS over-

683

estimate by about 0.05 the local AERONET measurements. 684 3.4.

3. Comparison with Maïdo measurements

	686
	687	the free troposphere are further compared to the measure-
	688	ments from the high altitude Maïdo Observatory (21.08 • S,
	689	55.38 • E; 2.2 km a.s.l). We use the measurements made
		by a Scanotron particle counter described in Foucart et al.

685

The in situ measurements on board the light plane in

Table A2

 A2 Same as TableA1but for F3 and F6. and origin of the submicron aerosol in the marine boundary layer 892 during the eastern mediterranean "sub-aero" experiment. Atmospheric 893Environment 40, 6245-6260. doi:10.1016/j.atmosenv.2006.03.059. Foucart, B., Sellegri, K., Tulet, P., Rose, C., Metzger, J.M., Picard, 895 D., 2018. High occurrence of new particle formation events at the 896 maïdo high-altitude observatory (2150 m), réunion (indian ocean). Gantt, B., Meskhidze, N., 2013. The physical and chemical characteristics 900 of marine primary organic aerosol: a review. Atmospheric chemistry Gao, R.S., Telg, H., McLaughlin, R.J., Ciciora, S.J., Watts, L.A., Richard-903 son, M.S., Schwarz, J.P., Perring, A.E., Thornberry, T.D., Rollins, A.W., Hamill, P., Giordano, M., Ward, C., Giles, D., Holben, B., 2016. An 915 aeronet-based aerosol classification using the mahalanobis distance. At-916 mospheric Environment 140, 213-233. doi:10.1016/j.atmosenv.2016.Holben, B., Eck, T., Slutsker, I., Tanré, D., Buis, J., Setzer, A., Vermote, 919 E., Reagan, J., Kaufman, Y., Nakajima, T., Lavenu, F., Jankowiak, I., Clerbaux, C., Coheur, P., Engelen, R.J., Errera, Q., Flemming, 925 J., George, M., Granier, C., Hadji-Lazaro, J., Huijnen, V., Hurtmans, 926 D., Jones, L., Kaiser, J.W., Kapsomenakis, J., Lefever, K., Leitão, J., 927 Razinger, M., Richter, A., Schultz, M.G., Simmons, A.J., Suttie, M., 928 Stein, O., Thépaut, J.N., Thouret, V., Vrekoussis, M., Zerefos, C., 929 2013. The MACC reanalysis: an 8 yr data set of atmospheric com-930 position. Atmospheric Chemistry and Physics 13, 4073-4109. doi:10. Karol, Y., Tanré, D., Goloub, P., Vervaerde, C., Balois, J.Y., Blarel, L., 933 Podvin, T., Mortier, A., Chaikovsky, A., 2013. Airborne sun photome-934 ter plasma: concept, measurements, comparison of aerosol extinction 935 vertical profile with lidar. Atmospheric Measurement Techniques 6. Dauhut, T., Ducrocq, V., Filippi, J.B., Gazen, D., Geoffroy, O., Gheusi, 949 F., Honnert, R., Lafore, J.P., Lebeaupin Brossier, C., Libois, Q., Lunet, 950 T., Mari, C., Maric, T., Mascart, P., Mogé, M., Molinié, G., Nuissier, 951O., Pantillon, F., Peyrillé, P., Pergaud, J., Perraud, E., Pianezze, J., Re-952 delsperger, J.L., Ricard, D., Richard, E., Riette, S., Rodier, Q., Schoet-953 ter, R., Seyfried, L., Stein, J., Suhre, K., Taufour, M., Thouron, O., 954 Turner, S.,Verrelle, A., Vié, B., Visentin, F., Vionnet, V., Wautelet,

				F3 (Group I) and F6
			FLEXPART simulation. Flight: F3
	897	Altitude 0-500 m	500-1500 m	> 1500 m Total	931
	898	Altitude release: 200 m Atmospheric Chemistry and Physics 18, 9243-9261. doi:10.5194/	5194/acp-13-4073-2013.
		Ocean acp-18-9243-2018. 43%	38%	15%	96%
		Island	2%	1%	1%	4%
	901		Altitude release: 1000 m
		Ocean and physics 13, 3979-3996. 19 % 42%	27%	88%	936
						doi:10.5194/amt-6-2383-2013.
		Island	2%	7%	4%	12%
			Altitude release: 1500 m
		Ocean	8%	41%	39%	88%
		Island	2%	8%	3%	13%
		FLEXPART simulation. Flight: F6 12h-simulation
			Altitude release: 1500 m
		Ocean	0.0%	35%	65%	100%
		Island	0.0%	0.0%	0.0%	0.0%
	917
		06.002.
	920
		Smirnov, A., 1998. Aeronet-a federated instrument network and data

sition 894 899 902 904 Markovic, M.Z., Bates, T.S., Johnson, J.E., Fahey, D.W., 2016. A light-905 weight, high-sensitivity particle spectrometer for pm2.5 aerosol mea-906 surements. Aerosol Science and Technology 50, 88-99. doi:10.1080/ 907 02786826.2015.1131809. 908 Guilpart, E., Vimeux, F., Evan, S., Brioude, J., Metzger, J.M., Barthe, C., 909 Risi, C., Cattani, O., 2017. The isotopic composition of near-surface wa-910 ter vapor at the maïdo observatory (reunion island, southwestern indian 911 ocean) documents the controls of the humidity of the subtropical tro-912 posphere. Journal of Geophysical Research: Atmospheres 122, 9628-913 9650. doi:10.1002/2017JD026791. 914 918 921 archive for aerosol characterization. Remote Sensing of Environment 922 66, 1-16. doi:10.1016/S0034-4257(98)00031-5. 923 Inness, A., Baier, F., Benedetti, A., Bouarar, I., Chabrillat, S., Clark, 924 H., 932 937 Koren, I., Dagan, G., Altaratz, O., 2014. From aerosol-limited to invigora-938 tion of warm convective clouds. Science 344, 1143-1146. doi:10.1126/ 939 science.1252595. 940 Kusmierczyk-Michulec, J., de Leeuw, G., Gonzalez, C.R., 2002. Empir-941 ical relationships between aerosol mass concentrations and Ångström 942 parameter. Geophysical Research Letters 29, 49-1-49-4. doi:10.1029/ 943 2001GL014128. 944 Lac, C., Chaboureau, J.P., Masson, V., Pinty, J.P., Tulet, P., Escobar, J., 945 Leriche, M., Barthe, C., Aouizerats, B., Augros, C., Aumond, P., Au-946 guste, F., Bechtold, P., Berthet, S., Bielli, S., Bosseur, F., Caumont, 947 O., Cohard, J.M., Colin, J., Couvreux, F., Cuxart, J., Delautier, G., 948 955 P., 2018. Overview of the meso-nh model version 5.4 and its applica-956 tions.

Geoscientific Model Development 11, 1929-1969. doi:10.5194/ 957 gmd-11-1929-2018.

In this paper, altitudes are given above sea level (a.s.l)

718

Overall, the MAGIC and TSI measurements made dur-

Acknowledgements

This work was supported by the French national pro-823 gramme LEFE/INSU. We are grateful to the Labex CAPPA 824 (ANR-11-LABX-0005-01) which has funded this work in 825 the context of the Cloud-aerosol interactions work package. ECMWF for providing freely reanaly-

Table 5 Comparison of average nighttime concentrations (Maïdo Observatory) and concentrations in the free troposphere during flights.

Mean nighttime concentration

Concentrations during flights Mean nighttime concentration at Maïdo D-1 (21∶00-03∶00 UTC) (in free troposphere) at Maïdo D+1 (21∶00-03∶00 UTC)

Scanotron TSI POPS TSI MAGIC Scanotron TSI Date (#/cm 3 ) (#/cm 3 ) (#/cm 3 ) (#/cm 3 ) (#/cm 3 ) (#/cm 3 ) (#/cm


Fabrice Boust
email: fabrice.boust@onera.fr

Thomas Lepetit

Shah Nawaz Burokur

Metagrating absorber: Design and implementation

The extent to which the introduction of subwavelength spatial modulation of electromagnetic properties improves absorption performances is studied. The proposed absorber represents an evolution from the Salisbury screen whereby the uniform resistive layer is replaced by a metagrating. A periodic supercell that supports only the specular reflection is first designed and load impedances are then engineered to suppress this diffraction mode. To experimentally demonstrate the concept, four prototypes are fabricated and tested in the microwave domain around 10 GHz. Furthermore, the performances assessed by a merit factor derived from Rozanov's bound show that the use of metagratings opens up good perspectives for improving the state of art. Our findings can pave the way towards the development of high performance absorbers for applications across a broad frequency spectrum. © 2020 Optica Publishing Group Electromagnetic (EM) absorbers are of prime importance for applications related to radar cross section reduction, EM interference, EM compatibility, EM protection [1] as well as antennas for mutual coupling reduction [2]. Absorbers are traditionally designed by combining homogeneous layers. The intrinsic characteristics of these layers, such as thickness, electric permittivity or magnetic permeability, are selected to best suit the intended application. Conventional absorbers include Dällenbach structures [3], Salisbury screens [4] and Jaumann layers [5][6]. A homogeneous magnetodielectric layer placed over a ground plane is used to dissipate incident power in Dällenbach absorbers. The Salisbury screen, composed of a resistive sheet with tuned resistance placed at a distance of a quarter wavelength over a ground plane, is used to generate destructive interference in a narrow frequency band. To widen the absorption bandwidth, several Salisbury screens can be stacked over each other to form the Jaumann layer. A bound on the thickness to bandwidth ratio of absorbers was established in 2000 by K. Rozanov [7]. It is important to expess the bandwidth in terms of wavelength since the major difficulty is to achieve absorption at low frequencies with a minimum thickness.

Today, absorbing structure designs are mostly oriented toward the use of frequency selective surfaces (FSS).1 As such, high impedance surfaces (HIS) have been added to thin Salisbury screens to minimize reflectivity over a wide frequency band [8][9] Metasurfaces where the imaginary parts of the electromagnetic parameters are manipulated to enhance the tangent losses have also been employed to increase losses to create an efficient absorber [10][11]. In these approaches, structuring is carried out on a small scale with respect to wavelength and with a single type of pattern to consider the layers as homogeneous and assign them a uniform permittivity or effective impedance. Under these conditions, the performance of these absorbers is still subjected to Rozanov's bound. To enhance the performance of absorbers, the circuit analog (CA) absorber, further characterized by reactive components in addition to resistive sheets, provides an effective solution to overcome the problems of narrow bandwidth and high thickness [12][13]. As such, the CA absorber exploits the use of FSSs or metasurfaces loaded with lumped resistors.

The concept of metagratings has evolved from classical diffraction gratings and represents an interesting alternative to metasurfaces for wavefront manipulation [14]. However, there is only a limited number of meta-atoms in a metagrating supercell compared to the densely packed metasurface supercell incorporating numerous meta-atoms with subwavelength arrangement. Therefore, unlike metasurfaces, the sparse arrangement of the constituting scatterers does not make possible a description of metagratings in terms of surface impedances. A metagrating, an array of scatterers separated by a distance of the order of the operating wavelength, enables diffraction engineering by canceling a number of undesired propagating diffraction orders and keeping only desired ones for radiation. As such, perfect anomalous reflection and perfect beam splitting [14][15][16][17][18][19][20] mode conversion [21], orbital angular momentum beams [22] and nearfield focusing [23][24] have been demonstrated with metagratings.

In this Letter, we propose to examine, to which extent, the modulation of electromagnetic properties at a scale smaller than the wavelength allows an increase in the absorption bandwidth.

Consequently, structures derived from the Salisbury screen by replacing the uniform resistive layer by a metagrating are studied here. The theoretical design is presented and validated by numerical simulations and experimental measurements for a 10 GHz operating frequency. The design of such metagrating absorber relies on the use of a periodic supercell that supports a single diffracted order, namely, the specular reflection and the aim is then to suppress this diffraction mode. The performances obtained from the proposed absorbing structure are further compared to a Salisbury screen of reference with the help of a merit factor derived from Rozanov's bound. A related study actually focuses on an analytical design methodology so as to show the influence of impedance density composition on the absorption bandwidth [25]. Such design concept is the key for future development of thin EM absorbers in various applications including, communication systems, stealth, camouflaging and noninvasive probing.

In a general manner, a metagrating is represented by a periodic array of supercells. Each supercell is itself composed of N thin wires. When the metagrating is illuminated by an incident plane wave, polarization line currents are excited in the wires, resulting in a scattered field represented by Floquet-Bloch modes defined by the period L of the supercell. As such, the diffraction angles of propagating orders m can be found via the grating formula: L(sinm -sini) = m, where m represents the number of an order and i is the incidence angle of the illuminating plane wave. A line current is represented by the 2-D Dirac delta function (y,z) that allows analytically calculating the scattered field composed of both propagating and evanescent diffraction orders.

As presented in Fig. 1, we exploit a reflective metagrating where the load wires are on the top face of a grounded dielectric substrate. The wires are excited by a TE-polarized plane wave illumination, where exp(it) is assumed, and as detailed in Ref. [18], the complex amplitude 𝐴 of the electric field of the m th diffracted plane wave is given by:

𝑨 𝒎 𝑻𝑬 𝒌 𝒌𝜼 𝟐𝑳 𝟏 𝑹 𝒎 𝑻𝑬 𝒆 𝒊𝜷 𝒎 𝒉 𝜷 𝒎 ∑ 𝑰 𝒒 𝒆 𝒊𝝃 𝒎 𝒒 𝟏 𝒅 𝜹 𝒎𝟎 𝑹 𝟎 𝑻𝑬 𝒆 𝟐𝒊𝜷 𝟎 𝒉 𝑵 𝒒 𝟏 (1)
where k and  are, respectively, the wavenumber and the characteristic impedance in air. m and m are the tangential and normal components of the plane waves wavevector, respectively, and are calculated as 𝜉 𝑘 sin 𝜃 2𝜋𝑚/𝐿 and 𝛽 𝑘 𝜉 . 𝑅 is the corresponding Fresnel's reflection coefficient and is given by:

𝑹 𝒎 𝑻𝑬 𝒊𝜸 𝒎 𝑻𝑬 𝐭𝐚𝐧 𝜷 𝒎 𝒔 𝒉 𝟏 𝒊𝜸 𝒎 𝑻𝑬 𝐭𝐚𝐧 𝜷 𝒎 𝒔 𝒉 𝟏 (2)
with 𝛾 𝜇 and 𝛽 𝜀 𝜇 𝑘 𝜉 . N = M line currents Iq must be exploited in a supercell to control the complex amplitudes of the M propagating diffraction orders in a metagrating. However, active and/or lossy impedances are required for a perfect control [19]. As such, if there are enough wires per period and if the line currents Iq are engineered, it is possible to control the electric field reflected by the metagrating. Iq can be tailored by engineering the load-impedance density Zq of each wire. These two quantities are related by Ohm's law, which is written as [26]:

𝒁 𝒒 𝑰 𝒒 𝑬 𝒙 𝒆𝒙𝒄 𝒚 𝒒 , 𝒉 𝒁 𝒊𝒏 𝑰 𝒒 ∑ 𝒁 𝒒𝒑 𝒎 𝑰 𝒑 𝑵 𝒑 𝟏 (3)
𝐸 𝑦 , ℎ 𝑒 𝑅 𝑒 𝑒 is the external electric field created by the incident wave and reflected from the grounded substrate and 𝑍 𝑘𝜂𝐻 𝑘𝑟 /4 is the input-impedance density of the wire. 𝑍 are the mutual-impedance densities corresponding to the interaction with the other wires and depend only on the geometry of the wires. The load-impedance density Zq can be calculated using the electromagnetic simulation-based local periodic approximation (LPA) approach proposed in Ref. [27]. Full-wave simulations are first performed on an elementary cell placed in an infinite array by applying periodic boundary conditions and illuminated by an incident plane wave to calculate the reflection coefficient. In such configuration, the cell size being of subwavelength scale limits reflection to the specular term. Assuming small diameter wires, the reflection coefficient obtained from simulations is then used to analytically calculate Zq, by subtracting the mutual impedance density due to effects of the interactions between wires and geometrical parameters such as substrate's thickness. It is thus possible to associate a three-dimensional design of an elementary cell and a load impedance, to be used in any metagrating.

In this study, we consider a Salisbury-like screen composed of a 3 mm thick grounded dielectric substrate of relative permittivity s = 2.2 and a metagrating layer with 1 line current per supercell in a first step, required to cancel the specular reflection at a single frequency. For the sake of simplicity, we consider normal incidence illumination. The length L of the supercell is chosen to be L <  such that only the 0 th diffracted order, which is the specular reflection, is present. Equation (1) can then be written as:

𝑨 𝟎 𝑻𝑬 𝒌 𝜼 𝟐𝑳 𝟏 𝑹 𝟎 𝑻𝑬 𝒆 𝒊𝒌𝒉 ∑ 𝑰 𝒒 𝑹 𝟎 𝑻𝑬 𝒆 𝟐𝒊𝒌𝒉 𝑵 𝒒 𝟏 (4)
with 𝑘 .

Solving analytically (3) allows expressing the current Iq as a function of Zq. The currents Iq is then injected into (4), which expresses the amplitude of the specular reflection in terms of the impedance densities Zq. Besides, our main goal is to further expand the absorption frequency bandwidth. Hence, several frequency points are targeted and a supplemental wire is engineered for each additional frequency by considering the limited number N of wires that can be placed in the same supercell of length L. Since (4) must be written for several frequencies, it is necessary to have an impedance model with frequency-independent parameters. Here, we consider a load impedance composed of a resistance and a capacitance:

𝒁 𝒒 𝒓 𝒒 𝒊 𝒄 𝒒 𝟏𝟎 𝟏𝟓 𝝎 (5)
where rq and cq are expressed in /m and fF.m, respectively. For N = 1, r1 et c1 can be found algebraically. For N = 2, 3 and 4, Eq. ( 4) is written for a set of frequencies 𝑓 chosen to maximize the bandwidth around 10 GHz with a magnitude of -15 dB. rq and cq are determined to reach a minimum for the quantity ∑ |𝐴 𝑓 | . Table I shows the values of rq and cq calculated for four structures, where each structure only differs by the number of wires N. The size L of the supercell is calculated to be 20 mm and the distance between two wires is L/N.

To engineer the calculated impedances, the elementary cell composed of two microstrip capacitances in series with a lumped resistance is selected, as illustrated in the photograph of the N = 1 fabricated metagrating sample in Fig. 2. The impedance densities are calculated for a set of values of resistance r between 10 and 100  and width A between 1 mm and 5 mm for two values of  (8 mm and 10 mm). The gap in the microstrip capacitances is fixed to 0.25 mm. By interpolation, the different widths and resistance values leading to the required impedances are determined. For the resistances, commercially available components are used. In the analytical calculations, the resistors are assumed to be perfect and the metal pads to be infinitely conductive, while in full-wave numerical simulations, a RLC model is used. A simple model composed of an inductance in series with the resistance and a capacitance in parallel is considered [28]. Table I summarizes the different values of  and A as well as the different components values r, Lr and Cr in the resistances. The period  along the xdirection is chosen so that the commercially available resistance values are not too far from the calculated ones, with a difference of less than 15%.

To verify the performances of the designed metagrating absorbers, full-wave simulations based on the finite-element method (FEM) technique using COMSOL Multiphysics are conducted. For experimental tests, reflection measurements are performed on samples having dimensions 140 mm × 140 mm in a microwave anechoic chamber dedicated to radar-cross-section bistatic measurements, where transmitting and receiving horn antennas are mounted on a circular track of 5 m radius. The transmitter is fixed at -2° and the receiver is fixed at +2°, as schematically illustrated in the inset of Fig. 2. The different reflection coefficients are shown in Fig. 3. To evaluate the obtained performances, we use a factor of merit derived from Rozanov's bound for non-magnetic materials. In Ref. [7], K Rozanov establishes that the thickness of a non-magnetic broadband absorbing screen with a reflection coefficient less than ΓdB on a continuous wavelength Δ is necessarily greater than:

𝒅 𝑹𝒐𝒛𝒂𝒏𝒐𝒗 𝐥𝐧 𝟏𝟎 𝟒𝟎𝝅 𝟐 𝚫𝝀 |𝚪 𝒅𝑩 | (6)
Consequently, the performance of the designed absorbing screens is evaluated by a factor of merit defined as the ratio of the Rozanov's thickness to the effective thickness of the screen:

𝑩𝑭𝑴 𝚪 𝒅𝑩 𝒅 𝑹𝒐𝒛𝒂𝒏𝒐𝒗 𝒅 𝐥𝐧 𝟏𝟎 𝟒𝟎𝝅 𝟐 |𝚪 𝒅𝑩 | 𝚫𝝀 𝒅 (7)
Such factor of merit, necessarily lower than unity for non-magnetic structures, makes possible a comparison of the performances of absorbing screens that operate at different frequency bands and for different levels of absorption, as shown in Table I.

As illustrated in Table I, the BFM is extracted for two different levels of absorption. The best performance is obtained for the experimental 3-wires metagrating absorber with a factor of merit of 0.26 compared to 0.16 for the analytically calculated Salisbury screen of similar thickness. A qualitative agreement is obtained between simulations and experiments provided that the resistance parasitic terms are taken into account in simulations. The slight residual frequency shift between 3D numerical simulations with parasitic terms and experiments is due to the model used for the resistors that does not take into account the welding of the pads and also to the dispersion of the components. In addition, the LPA method used to calculate impedance densities takes into account the interactions between wires assuming that they are equivalent to infinitely thin wires. This approximation is less and less valid as the number of wires in a supercell increases. This is one of the hypotheses with the dispersion of resistance characteristics that could explain the poor performance of the 4-wires absorber. In conclusion, the impedance modulation induced by the designed metagrating structures achieves a better performance than the Salisbury screen, particularly when high level of absorption is desired. For a relevant comparison with widely developed metasurface-based absorbers, improvements on the proposed metagrating are necessary, such as by stacking several layers and by considering an elaborated load impedance model. In this work, the proposed metagrating abosrber concept is validated for a single polarization, but can be further extended to two simultaneous polarizations by combining two perpendicular wire arrays. Finally, the approach can be extended to metagratings with several diffracted orders and also to multilayer systems, in analogy to a Salisbury screen with a front-facing impedance matching layer, which would possibly lead to wider absorption bandwidth. It would be then necessary to replace the lumped resistors by resistive layers to facilitate their integration. Disclosures. The authors declare no conflicts of interest.

Fig. 1 .

 1 Fig. 1. Schematic principle of the metagrating absorber implemented by a single load wire per period (N = 1).

Fig. 2 .

 2 Fig. 2. Photograph of a fabricated sample composed of one wire per supercell (N = 1). The geometrical dimensions are: L = 20 mm,  = 10 mm and A = 2.2 mm. The gap in the microstrip capacitances is 0.25 mm. The resistance value is 50 .

Fig. 3 .

 3 Fig. 3. Performances in terms of reflection coefficients for different metagrating structures composed of N = 1, 2, 3 and 4 load wires. Theoretical calculations, numerical simulations and experimental measurements are shown and compared to the performance of a Salisbury screen of similar thickness.

Table I . Parameters and performances of metagrating structures with different number of wires.

 I

				Cell design						BFM (Analytical)	BFM (Experimental)
	N	rq (/m)	cq (fF.m)	(capacitances and resistances)	 (mm)	A (mm)	r ()	Lr (nH)	Cr (fF)	-10 dB	-15 dB	-10 dB	-15 dB
	1	8530	0.26	2 C + 1 R	10	2.2	50	0.9	100	0.03	0.024	0.052	0.027
	2	11190 8790	0.20 0.29	2 C + 1 R 2 C + 1 R	8 8	2.4 3.5	50 50	0.9 0.9	100 100	0.13	0.15	0.15	0.15
		3860	0.17	2 C + 1 R	8	1.8	18	0.4	60
	3	10420	0.32	2 C + 1 R	8	3.85	50	0.9	100	0.15	0.18	0.21	0.26
		12210	0.22	2 C + 1 R	8	2.6	50	0.9	100
		10885	0.33	2 C + 1 R	8	4.0	50	0.9	100
	4	5340 23380	0.14 0.21	2 C + 1 R 2 C + 1 R	8 8	1.5 2.7	25 100	0.45 0.9	70 100	0.25	0.33	0.12	0
		120	0.14	1 C	8	4.2	-	-	-
	Salisbury									0.20	0.16

Data availability. Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.


Alesia Herasimenka
email: alesia.herasimenka@univ-cotedazur.fr.

Lamberto Dell'elce

Jean-Baptiste Caillau

Jean-Baptiste Pomet

Controllability Properties of Solar Sails

Trajectory design and station keeping of solar sails about a celestial body can be formulated as control problems with positivity constraints. Specifically, when re-emitted radiation is neglected and the sail is modeled as a flat surface, which are reasonable assumptions for control purposes, force generated by solar radiation pressure is contained in a pointed convex cone of revolution with axis towards the Sun-satellite direction.

Therefore, classical approaches to infer controllability based on the Lie algebra rank condition do not apply to these problems. This study offers a novel condition to decide on controllability of control systems with positivity constraints. This condition is effective as it can be verified by solving an auxiliary convex optimization problem for which reliable numerical methods are available. A crucial ingredient of this approach is the theory of positive trigonometric polynomials. The practical interest of this condition is the assessment of a minimum requirement on the optical properties of the sail, which may be of use for mission design purposes.

I. Introduction

Solar sails offer a propellant-less solution to achieve interplanetary transfers, planet escapes, and de-orbiting maneuvers by leveraging on solar radiation pressure (SRP) [START_REF] Mcinnes | Solar Sailing[END_REF]. Although very few solar sail missions were launched, the possibility to use SRP as an inexhaustible source of propulsion attracted the interest of researchers since decades, and several contributions on the guidance and control of solar sails are available. Specifically, a large body of literature focuses on the mathematical formulation and numerical solution of optimal interplanetary transfers (mostly minimum time) using optimization techniques [START_REF] Macdonald | Heliocentric Solar Sail Orbit Transfers with Locally Optimal Control Laws[END_REF], indirect methods [START_REF] Gong | Solar sail time-optimal interplanetary transfer trajectory design[END_REF] or even neural networks [START_REF] Dachwald | Optimal Solar Sail Trajectories for Missions to the Outer Solar System[END_REF].

In addition, several contributions investigate locally optimal maneuvers, i.e., maximization of the instantaneous rate of change of a desired orbital element, in particular increasing semi-major axis for orbit-raising [START_REF] Mengali | Near-Optimal Solar-Sail Orbit-Raising from Low Earth Orbit[END_REF] or decreasing perigee altitude for de-orbiting applications [START_REF] Colombo | Modulating Solar Sail Control for End-Of-life Disposal with Solar Sails[END_REF]. Classical feedback algorithms are also used to find sub-optimal trajectories, as for instance Q-law in [START_REF] Niccolai | Solar sail heliocentric transfers with a Q-law[END_REF]. Direct methods are often preferred to tackle the numerical solution of optimal control problems (OCP) of solar sail transfers [START_REF] Caruso | Optimal solar sail trajectory approximation with finite Fourier series[END_REF] because they do not require an initial guess of adjoint variables, differently of indirect techniques, which were used in very few studies like [START_REF] Kim | Symmetries in the Optimal Control of Solar Sail Spacecraft[END_REF] and [START_REF] Sullo | Low-Thrust to Solar-Sail Trajectories: A Homotopic Approach[END_REF]. Surprisingly enough, very few studies on the controllability of solar sails are available to date, although an analysis of the reachable set of passive sails was discussed in [START_REF] Pellegrino | Reachability of a Passive Solar Sail in Earth Orbit[END_REF]. Most often, solutions of two-point boundary value problems coming for instance from optimal control are investigating without assessing beforehand whether the targeted point is within the reachable set of the control system from the initial point; exhibiting a solution of the two-point boundary value problem obviously proves reachability as a side result, but it is better to be able to certify non controllability and seek for solutions only when controllability holds.

A major difficulty in assessing the controllability of an SRP-actuated system is that the sail cannot generate a force with a positive component toward the direction of the Sun, so that classical tools of geometric control theory cannot be used. For example, Lie algebra of the system is full rank (unless a fully absorptive model of the sail is considered), but although it should indicate that the system is weakly controllable, theory requires that the control set is a neighborhood of the origin, i.e., both positive and negative controls should be generated, so that it is not sufficient to analyze the sailing problem. This aspect is particularly critical when considering a non-ideal sail (by ideal, we mean a perfectly reflective flat surface) because the control set is contained inside a strictly convex cone of revolution, whose angle depends on the optical properties of the sail.

The main contribution of the paper consists in a controllability check for non-ideal solar sails in a planet-centered orbit. This requirement is aimed at assessing whether the sail at hand is capable of decreasing or increasing all possible functions of the Keplerian integrals of motion over an orbital period. Given some optical properties, a convex cone containing all possible directions of the SRP is first defined. Then, the condition is verified by means of a worst-case optimization problem characterized by a finite number of design variables and a two-parameter family of inequality constraints, namely, the clock angle of the convex cone and the true anomaly of the sail. Numerical solution of this semi-infinite problem is achieved by leveraging on the formalism of squared functional systems [START_REF] Nesterov | Squared Functional Systems and Optimization Problems[END_REF][START_REF] Dumitrescu | Positive Trigonometric Polynomials and Signal Processing Applications[END_REF] to exactly enforce inequality constraints for all values of the true anomaly and the clock angle. No discretization is done to solve the problem numerically. Eventually, the semi-infinite problem is recast into a finite-dimensional convex programming with a finite number of linear matrix inequalities (LMI) and an unique well-defined solution. Non-satisfaction of the condition entails some local non-controllability of the system for the given value of the cone angle (and, consequently, of the optical properties) and orbital conditions. Hence, a fine analysis covering (as much as possible numerically) the entire phase space of orbital elements is carried out to determine the minimum cone angle for a large range of orbits. It is shown that a universal (namely, planet independent) minimum angle exists that satisfies the condition for all orbits. Its value is about 60 degrees (note that 0 and 90 degrees correspond to fully absorptive and perfectly reflective sails, respectively). The result indicates that the sail does not have to be ideal to satisfy the requirement.

The methodology is unaware of the specific source of non-ideality of the sail (e.g., specular or diffuse reflection, re-emitted radiation [START_REF] Mcinnes | Solar Sailing[END_REF]) since it only uses the conical hull of the control set. This result can be used to provide insight into the controllability of the sail during its lifetime, owing to the degradation of its optical properties discussed in [START_REF] Dachwald | Impact of Optical Degradation on Solar Sail Mission Performance[END_REF], and may support the design of real-life missions by serving as a minimal requirement to be satisfied.

Section II introduces the dynamical model of solar sails, geometry of the problem, equations of motion, and various assumptions used in this work. Section III outlines the novel condition for investigating local controllability of systems with peculiar constraints on the control set that apply to solar sails. An efficient numerical methodology based on convex programming to evaluate the aforementioned condition is detailed in section IV. Finally, this methodology is extensively used in Section V to deduce minimal requirements on the optical properties of solar sails. Controllability of heliocentric orbits is briefly discussed as well.

II. Dynamics of a solar sail

A. Force model Solar sails use SRP as propulsive means to control their trajectory. SRP is due to the interaction between photons and surface of the sail. The magnitude of the pressure depends on the Sun-sail distance, r (in AU). Specifically, denoting by Φ SR ≈ 1367 W m -2 the solar flux at r ⊕ = 1 AU and by c the speed of light, a simple model is [START_REF] Montenbruck | Satellite Orbits[END_REF]Chap. 3]:

P SR = Φ SR c r ⊕ r 2
A flat sail with surface A and mass m is considered in this work. The resulting force depends on various optical and geometrical properties of the sail and it is obtained by summing up contributions of the incoming, reflected, and thermal radiations, namely f a , f r , and f e . In addition, the reflected force is divided into specular and diffuse components, f rs and f ru , respectively. The first one is caused by photons that are reflected symmetrically with respect to the normal of the sail and create moment in the opposite direction. Conversely, diffuse reflection stems from surface roughness, which causes photons to be uniformly reflected in all directions, yielding a component of the force toward the vector normal to the sail. Finally, absorbed photons are then re-radiated in all directions with energy dependent on the temperature of the sail, generating another component of the force that is orthogonal to its surface. Figure 1a shows the directions of the various components. Denoting by ŝ the direction of the Sun, n the unit vector normal to the sail with positive projection toward ŝ, i.e., cos β := n • ŝ ≥ 0, and

t := n × ŝ n × ŝ × n = ŝ -cos β n sin β
the tangent unit vector in the plane generated by ŝ and n, the components of the specific force are [START_REF] Rios-Reyes | Generalized Model for Solar Sails[END_REF] f a = ε(r ) cos β(cos β n + sin β t)

f rs = ε(r ) ρ s cos β(cos β n -sin β t) f ru = ε(r ) B f ρ(1 -s) cos β n f e = ε(r ) (1 -ρ) ε f B f -ε b B b ε b + ε f cos β n (1)
(a) Schematic representation.

f SRP = f a + f rs + f ru + f e .

B. Parametrization of the control set

Controlling the sail attitude, i.e. the normal vector n, allows to change the direction and magnitude of the resulting SRP. A reliable inference of optical coefficients is indeed mandatory to accurately estimate the mapping between n and f SRP .

To carry out controllability analysis, solar sail dynamics is conveniently modeled as a nonlinear control-affine system (see Section II C), where the control variable is homogeneous to the force vector, namely u :=

f SRP ε(r )
. Control set U ⊂ R 3 is then given by:

U = f SRP (n) ε(r ) , n ∈ R 3 , n = 1 .
Figure 1b shows the projection of U on the plane generated by n and ŝ for various optical properties. The set is a surface of revolution with axis ŝ, and it is non-convex unless ρ = s = 1.

Note that the interior of the surface is not part of U . When re-emitted radiation is neglected, which is most often a reasonable assumption for control purposes, U contains the origin but mapping between n and u is non-smooth at this point. Two extreme cases can be identified: ideal sails are constituted by perfectly reflective surfaces (ρ = s = 1), whereas perfectly-absorptive surfaces are the worst-case scenario (ρ = 0, f e neglected) because SRP is systematically parallel to ŝ. Although sails are designed to be as ideal as possible, partial absorption of the energy is unavoidable in real-life applications and, in addition, optical properties exhibit degradation with time. Hence, the fraction of reflected radiation decreases with lifetime of the satellite, as discussed in [START_REF] Dachwald | Impact of Optical Degradation on Solar Sail Mission Performance[END_REF].

C. Equations of motion

Controllability of a solar sail in orbit around a celestial body is considered in this study. The following assumptions are introduced:

1. Orbital period of the sail is much smaller than the one of the heliocentric orbit of the attractor, so that variations of the Sun direction ŝ over a single orbit of the sail are neglected.

2. Solar eclipses are neglected. Targeting a certification of non-controllability in Section III, this assumption is conservative, and it has the major advantage of yielding results that are independent of the semi-major axis of the orbit.

3. Re-emitted radiation is neglected. In fact, this component of SRP can be reasonably regarded as a disturbance for control purposes.

Equations of motion are written in a set of Keplerian-like orbital elements, which leverages on the axial symmetry of the problem with respect to the Sun's direction. Namely, consider a reference frame with origin at the center of the planet, X axis towards ŝ, whereas Ŷ and Ẑ are arbitrarily chosen and form a right-hand frame. Because this study focuses on short-time controllability (time is of the order of one orbital period), motion of this frame is neglected by virtue of the first assumption above. Let γ 1 , γ 2 , γ 3 be Euler angles orienting the eccentricity vector according to a X -Y -X rotation as depicted in Fig. 2, and a, e, and f be semi-major axis, eccentricity and true anomaly, respectively. Motion of slow elements,

I = (γ 1 , γ 2 , γ 3 , a, e) T ∈ M , is governed by d I d t = ε(r ) a 1 -e 2 µ G(I, f ) R(I, f ) u d f d t = ω(I, f ) + ε(r )F (I, f ) R(I, f ) u (2)
where components of u are in the reference frame, R(

I, f ) = R X (γ 3 + f )R Y (γ 2 )R X (γ 1
) is the rotation matrix from reference to local-vertical local-horizontal frames 5 , ω(I, f ) = µ a (1 -e 2 ) 3 (1+ e cos f ) 2 , both F (I, f ) and G(I, f ) can be deduced from Gauss variational equations (GVE) of classical elements, where G(I, f ) is:

G =                  0 0 sin (γ 3 + f ) sin γ 2 (1 + e cos f ) 0 0 cos (γ 3 + f ) 1 + e cos f - cos f 7 te 2 + e cos f 1 + e cos f sin f e cos (γ 3 + f ) 1 + e cos f 2 a e 1 -e 2 sin f 2 a e 1 -e 2 (1 + e cos f ) 0 sin f e cos 2 f + 2 cos f + e 1 + e cos f 0                 
The peculiar choice of Euler angles follows from the symmetry of System (2), and it has the main consequence that controllability results in Section V are independent of γ 1 . We also note that

(1 + e cos f ) G(I, f )R(I, f
) is a trigonometric polynomial in f . This has significant advantages for the numerical methodology detailed in Section IV.

Finally, orbital perturbations (other than SRP) are not included in Eq. ( 2) because we are interest in investigating geometric obstructions to the controllability of solar sails regardless their size.

III. Controllability of a solar sail

We are interested in studying obstructions to controllability of solar sails in orbit about a celestial body. More precisely, we want to assess the existence of variations of the current orbital elements set that a sail cannot generate after a single orbital period.

The classical approach to inspect controllability of control-affine systems with periodic drift 5 Here, R A (ϕ) denotes the rotation matrix of angle ϕ about the axis Â.

is detailed in [START_REF] Bonnard | Contrôlabilité des systèmes non linéaires[END_REF] and [START_REF] Jurdjevic | Geometric Control Theory[END_REF]Chap. 4]. Specifically, global controllability (that is existence of an admissible control steering the system from any initial point towards any target) is guaranteed (sufficient condition) provided:

1. the drift of the system is periodic (or, more generally, recurrent);

2. the Lie algebra rank condition (LARC) condition holds, namely the set of vector fields defining the control-affine system are bracket generating;

3. the convex hull of the control set U is a neighborhood of the origin in R m .

Keplerian motion satisfies the first condition. Two different cases are distinguished for the computation of the rank of Lie algebra: First, perfectly-absorptive sails (i.e., ρ = 0) are such that the control set degenerates to a segment aligned to ŝ, as shown in Fig. 1b. Appendix A offers a detailed evaluation of Lie brackets in this case, which shows that the algebra is rank deficient. This allows to conclude the non-controllability of the system and to deduce an integral of motion, namely the projection of the angular momentum towards ŝ [19, Chap. 12]. Second, real-life sails are such that ρ is strictly positive and smaller than 1. In this case, control set is no more degenerate and the system is bracket generating. Nevertheless, controllability of this system cannot be assessed because U does not contain the origin in its interior (the origin is on the boundary of the control set as depicted in Fig. 1b), so that the third condition above is not satisfied, and the classical approach cannot be used to investigate controllability of this system.

A novel sufficient condition relaxing the third requirement was proposed in [START_REF] Herasimenka | Controllability Test for Fast-Oscillating Systems with Constrained Control. Application to Solar Sailing[END_REF] and [START_REF] Caillau | On the Controllability of Nonlinear Systems with a Periodic Drift[END_REF]. System ( 2) is controllable wrt. to the slow variables, meaning that there exists an admissible control allowing to reach any final orbit from any initial one, 6 if the control set U contains the origin (even only on its boundary), and if for all I in M the following condition holds:

cone d I(f, I, u) d t , u ∈ U, f ∈ S 1 = T I M. (3)
Here before the operator cone indicates the conical hull: for a subset A of a vector space, cone(A)

is the set of linear combinations with nonnegative coefficients of vectors of A,

k i=1 λ i x i , k ∈ N, λ 1 , . . . , λ k ≥ 0.
It is always a convex pointed (containing the origin) cone. This result is proved in [START_REF] Herasimenka | Controllability Test for Fast-Oscillating Systems with Constrained Control. Application to Solar Sailing[END_REF][START_REF] Caillau | On the Controllability of Nonlinear Systems with a Periodic Drift[END_REF]. If

(3) only holds at one point, it remains true around it and controllability in the corresponding neighborhood follows. Conservely, (3) cannot hold if, for some orbit I 0 , there exists a nonzero one-form p I 0 ∈ T * I0 M such that

p I 0 , d I(f, I 0 , u) d t > 0, f ∈ S 1 , u ∈ U. (4)
In this case, there is an obstruction to local controllability in the following sense: one can find a neighborhood of I 0 on which motion is prescribed to a half-space (see Fig. 3); orbits in the forbidden neighborhood half might be reached (e.g. if global controllability holds), but only by trajectories that must leave this neighborhood. For a detailed discussion on the gap between (3) and (4), we refer to [START_REF] Caillau | On the Controllability of Nonlinear Systems with a Periodic Drift[END_REF].

IV. Numerical methodology to inspect local controllability

Given some optical properties of the sail and orbital state I, we are interested in determining if Eq. ( 4) has any non-trivial solution p I = 0. Two manipulations are introduced to facilitate this task.

not play any role; any initial/final longitude can be departed from/reached, up to some additional time spent along the periodic initial/final orbit. First, the time derivative of I in Eq. ( 4) is replaced by G(I, f )u, where

G(I, f ) := (1 + e cos f )G(I, f )R(I, f ),
which is a second-degree trigonometric polynomial in f . Such property offers major benefits when positivity constraints are numerically enforced in Section IV B. This operation has no impact on the sign of Eq. ( 4). We also note that System (2) is axially symmetric with respect to the Sun-planet direction, and that semi-major axis and planetary constant have no impact on the sign of Eq. ( 4).

Hence, all outcomes of this controllability study are independent of both the semi-major axis and γ 1 (because of symmetry), and they are valid for any attractor (spherical symmetric central body), since magnitude of SRP does not impact the non-controllability condition (which is a geometric obstruction).

Second, control set U is replaced by its conical hull, K α := cone(U ), which is a cone of revolution of angle α, as illustrated in Fig. 4. This approximation makes the problem convex (in fact, convex programming is used to get numerical certifies of the feasibility problem in Section IV A). Replacing U by K α has no impact on the closure of the reachable set of the control system, as discussed in [START_REF] Caillau | On the Controllability of Nonlinear Systems with a Periodic Drift[END_REF].

Therefore, non-controllability of the system with controls in K α implies non-controllability of the system with the original control set, U . Neglecting thermal radiation (this simplification is not strictly necessary), cone angle α can be directly deduced from the optical properties of the sail introduced in Eq. ( 1). The relation is obtained by solving:

tan α = max β∈[0, π 2 ] f SRP • ŝ I -ŝŝ T f SRP = max β∈[0, π 2 ] ρ s sin 2β + B f ρ(1 -s) sin β 1 + ρ s cos 2β + B f ρ(1 -s) cos β . (5)
This condition holds for:

β * = cos -1   -B f ρ(1 -s)(3ρs + 1) + B 2 f ρ 2 (1 -s) 2 ((3ρs -1) 2 -4ρs) -32ρ 2 s 2 (ρs -1) 8ρs   .
If B f = 0, Eq. ( 5) simplifies to

α(ρ, s) = tan -1 ρ s 1 -ρ 2 s 2 , ρ s = tan α √ 1 + tan 2 α • (6)
Hence, Eq. ( 4) is finally recast into if ∃ p I ∈ T * M, p I = 0 such that

p I , G(I, f ) u > 0, f ∈ S 1 , u ∈ K α . (7)
A. Constructive approach to verify the controllability condition A practical check of the feasibility Problem ( 7) is carried out by solving the auxiliary optimization problem max

J, p I ≤1 J s.t. p I , G(I, f ) u ≥ J, f ∈ S 1 , u ∈ ∂K α , u = 1. (8)
The constraint p I ≤ 1 is preferred to p I = 1 to preserve convexity of Problem [START_REF] Caruso | Optimal solar sail trajectory approximation with finite Fourier series[END_REF]. Problem ( 8) is convex and semi-infinite, because inequality constraints need to be enforced on two infinite sets, namely for all true anomalies between 0 and 2π and for all u on the surface of the cone. Evaluating inequalities in the interior of the cone is not necessary because dynamics is affine in u. If J * , solution of Problem [START_REF] Caruso | Optimal solar sail trajectory approximation with finite Fourier series[END_REF], is positive, ( 7) is verified: then, as discussed in the previous section, for the cone angle α, there is an obstruction to local controllability around the orbit I.

Conversely, when (3) holds at I7 , both J * and the associated minimizer p I must be zero.

To emphasize the practical interest of Problem (8), Figure 5 depicts an example solution for a given orbit as a function of α (the detailed numerical algorithm to achieve these solutions is provided in Section IV B). One can check that a minimum cone angle exists (for the specific I used in this simulation) for which [START_REF] Niccolai | Solar sail heliocentric transfers with a Q-law[END_REF] (and the corresponding obstruction to local controllability) cannot hold. This angle can be mapped into minimal requirements for the reflectivity of the sail via Eq. ( 6), and it can be evaluated by solving

min α α s.t. J * (α) = 0 (9)
where J * (α) denotes solution of Problem (8) for a given α.

B. Optimization problem

Numerical solution of Problem ( 8) is achieved by using the formalism of positive trigonometric polynomials [START_REF] Nesterov | Squared Functional Systems and Optimization Problems[END_REF][START_REF] Dumitrescu | Positive Trigonometric Polynomials and Signal Processing Applications[END_REF] to enforce positivity constraints for all values of f and u without introducing any relaxation or discretization of the problem.

Let δ be an angle parametrizing control vectors on the surface of the cone, as shown in Fig. 6, Fig. 6: Parametrization of the control vector.

namely

u =         cos α cos δ sin α sin δ sin α        
Positivity of the following constraints must be verified:

p T I G(I, f )u -J ≥ 0 f ∈ S 1 , u ∈ ∂K α . (10)
Inspection of G(I, f )u reveals that Eq. ( 10) is a bivariate trigonometric polynomial of second degree in f and first degree in δ. Let •, • H be the Hermitian product of two complex-valued vectors,

(f, δ) = 1, e iδ T ⊗ 1, e if , e 2if T =
1, e if , e 2if , e iδ , e if e iδ , e 2if e iδ T the basis of bivariate trigonometric polynomials of degree 2 in f and 1 in δ, respectively (here, ⊗ denotes Kronecker's product). The left-hand term of Eq.( 10) can be reformulated as

p I , G(I, f ) u -J = p T I 1 l=-1 2 k=-2 Gu (k,l) e i k f e i l δ -J = Φ(f, δ), Gu p I -e 1 J H
where Gu (k,l) (I) is the kl-th coefficient of the Fourier transform8 of G(I, f )u, and e 1 = [1, 0, 0, 0, 0, 0] T .

The formalism of squared functional systems outlined in [START_REF] Nesterov | Squared Functional Systems and Optimization Problems[END_REF] and [START_REF] Dumitrescu | Positive Trigonometric Polynomials and Signal Processing Applications[END_REF] allows to recast the continuous positivity constraints into LMI. The corresponding squared functional system of Φ(f, δ)

is S 2 (f, δ) = Φ(f, δ)Φ H (f, δ), where Φ H (f, δ) denotes conjugate transpose of Φ(f, δ).
Let N be the dimension of Φ(f, δ) (6 in our application) and Λ H : C N → C N ×N be a linear operator mapping coefficients of polynomials in Φ(f, δ) to the squared base, so that application of Λ H on Φ(f, δ) yields

Λ H (Φ(f, δ)) = Φ(f, δ)Φ H (f, δ)
and define its adjoint operator Λ * H :

C N ×N → C N as Y, Λ H ( Gu) H ≡ Λ * H (Y ), Gu H , Y ∈ C N ×N , Gu ∈ C N .
Theory of squared functional systems postulated by Nesterov [START_REF] Nesterov | Squared Functional Systems and Optimization Problems[END_REF] proves that trigonometric polynomial is non-negative if and only if a Hermitian positive semidefinite matrix Y exists such that Gu = Λ * H (Y ). Dumitrescu extends this theory for multivariate trigonometric polynomials in [START_REF] Dumitrescu | Positive Trigonometric Polynomials and Signal Processing Applications[END_REF]Chap. 3] and shows that all nonnegative bivariate trigonometric polynomials can be written as sum-of-squares. This equivalence is false for three or more variables.

Thus, Φ(f, δ), Gu H is non-negative for all f ∈ S 1 and for all u ∈ K α if and only if a Hermitian positive semidefinite matrix Y exists such that Gu = Λ * H (Y ), namely

Φ(f, δ), Gu H ≥ 0, f ∈ S 1 , u ∈ K α ⇐⇒ ∃Y 0 : Gu = Λ * H (Y ).
In fact, it holds in this case that

Φ(f, δ), c H = Φ(f, δ), Λ * H (Y ) H = Λ H (Φ(f, δ)), Y H , = Φ(f, δ)Φ H (f, δ), Y H = Φ H (f, δ) Y Φ(f, δ) ≥ 0.
For trigonometric polynomials Λ * is given by

Λ * H (Y ) =                 tr( Y, T 00 ) . . . tr( Y, T kl ) . . . tr( Y, T 21 )                 k = 0, 1, 2, l = 0, 1.
where T j j = 0, 1, 2 are the elementary Toeplitz matrices with ones on the j-th diagonal and zeros elsewhere and T kl are obtained from a Kronecker product of such matrices, e.g.,

T 0 =     1 0 0 1     , T 1 =         0 1 0 0 0 1 0 0 0         , T 10 = T 0 ⊗ T 1 =                     0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0                    
Finally, the inequality in Eq. ( 10) is rewritten as an LMI:

p I , G(I, f ) u -J ≥ 0, f ∈ S 1 , u ∈ ∂K α ⇐⇒ ∃ Y 0 such that Gu p I -e 1 J = Λ * H (Y )
where Y ∈ C 6×6 is a Hermitian matrix to be determined. Hence, the finite-dimensional counterpart of Problem ( 8) is

min J, p I ≤1, Y ∈C 6×6 J s.t.: Y 0 Λ * H (Y ) = Gu p I -e 1 J (11)
Eventually, solution of Problem ( 9) is carried out by means of a simple bisection algorithm. The CVX software [START_REF] Grant | CVX: Matlab Software for Disciplined Convex Programming, version 2.1[END_REF][START_REF] Grant | Graph implementations for nonsmooth convex programs[END_REF] is used to solve the convex Problem [START_REF] Pellegrino | Reachability of a Passive Solar Sail in Earth Orbit[END_REF]. Fourier coefficients of G(I, f ) are evaluated by means of the fast Fourier transform (FFT) algorithm.

We stress that there is no relaxation of Problem [START_REF] Pellegrino | Reachability of a Passive Solar Sail in Earth Orbit[END_REF] with respect to Problem [START_REF] Caruso | Optimal solar sail trajectory approximation with finite Fourier series[END_REF]. Remarkably, enforcement of the constraint for all values of f and u is exact and stems from the trigonometric nature of G(I, f )u. Hence, this methodology outperforms the algorithm that we proposed in [START_REF] Herasimenka | Controllability Test for Fast-Oscillating Systems with Constrained Control. Application to Solar Sailing[END_REF],

where controls necessary to move towards the vertexes of a simplex encompassing I where explicitly evaluated. In that case, a large number of harmonics was necessary to model the control variable and, eventually, to achieve a conservative estimate of the minimal cone angle α, as depicted in Fig. 7.

V. Minimal optical requirements

Figure 8 shows the minimum cone angle satisfying the condition as a function of γ 2 and γ 3 for various values of eccentricity (we recall that semi-major axis and γ 1 have no influence on this angle). The minimal angle is symmetric with respect to γ 2 = 90 deg because p I , G(e, γ 2 , γ 3 , f ) u = -p I , G(e, π -γ 2 , γ 3 , f ) u . Solution is independent of γ 3 for circular orbits, as expected. Sensitivity with respect to γ 3 remains moderate even for larger eccentricities. The minimal angle approaches zero as sin(γ 2 ) → 0. In this case, ŝ is aligned with the angular momentum of the orbit.

On the other hand, for γ 2 = 90 deg, the Sun is in the orbital plane. Fig. 7: Convergence of the results obtained in [START_REF] Herasimenka | Controllability Test for Fast-Oscillating Systems with Constrained Control. Application to Solar Sailing[END_REF] as the number of harmonics increases (in red)

towards the minimum angle obtained with the methodology detailed in this paper (in blue).

Figure 9 represents α min as function of γ 3 or γ 2 (Fig. 9a and 9b, respectively) for various values of eccentricity. Results confirm high dependency of α min on γ 2 , and γ 3 for large eccentricity. Hence, controllability of near circular orbits requires more reflective sails with respect to high-eccentric orbits. Finally, we stress that the minimum angle α exists for all orbits, and it is systematically smaller than 90 deg, which means that the sail has not to be ideal to make System (2) controllable.

To compare with a real solar sail, optical properties of the NASA reference model [START_REF] Heaton | An Update to the NASA Reference Solar Sail Thrust Model[END_REF] (designed to support NEA Scout and Lunar Flashlight solar sail missions) correspond to a cone angle of 58.6 deg.

This value is sufficient to satisfy the proposed condition for most planet-centered orbits, except for highly inclined ones.

A. A comment on heliocentric orbits

Consider now a sail in a heliocentric orbit. This scenario can be a case for interplanetary trans- First, the rotation matrix R in Eq. ( 2) is removed, since the local vertical local horizontal frame is used, and ŝ is aligned with the radial direction. Moreover, the problem has central symmetry, so that results do not depend on any orbital element except for the eccentricity. For a perfectly absorptive solar sail, the dynamical system is not bracket generating, because the control is radial, as proved in [START_REF] Bonnard | Geometric optimal control of elliptic Keplerian orbits[END_REF]. The integral of motion related to this rank deficiency is the magnitude of the angular momentum. For a non-ideal sail, the system becomes bracket generating as soon as a Using the methodology of Section IV indicates that even a very poorly reflective sail (i. e., 0 < ρ 1) is locally controllable over an orbital period. However, an orbital period of heliocentric orbits is of the order of the year, this is why the proposed methodology is not well suitable to tackle realistic mission scenario of interplanetary transfers using solar sails. In fact, a sail can be moved towards any direction of the tangent manifold of the orbital elements set, but the requirement is unaware of the orbital period's magnitude, so that these maneuvers could possibly last for years.

VI. Conclusion

We introduced a novel requirement to assess the controllability of solar sails in orbit around a celestial body. This was formulated as a geometric condition on the reachable set of orbits (slow variables of the system, in contrast with longitude which is the fast one). A failure of this requirement typically indicates that a half neighborhood of the state vector exists where the sail cannot be moved. A numerical algorithm based on convex optimization and leveraging on the theory vectors, respectively, System (A2) can be rewritten as:

ẋ = F 0 (x) + ε uF 1 (x)
with F 0 the recurrent drift and F 1 the SRP perturbation.

             d r d t = v d v d t = - µ r 3 r + ε(r ) ŝu (A2)
where r = r . System (A2) provides two vector fields:

F 0 = v X ∂ ∂r X + v Y ∂ ∂r Y + v Z ∂ ∂r Z - r X r 3 ∂ ∂v X - r Y r 3 ∂ ∂v Y - r Z r 3 ∂ ∂v Z F 1 = s X ∂ ∂v X + s Y ∂ ∂v Y + s Z ∂ ∂v Z
To simplify, let us denote vector fields: All subsequent iterations are linear combinations of the previous vector fields. Thus, Lie algebra of the system (A2) has 5 independent vector fields if ŝ • r = 0:

v ∂ ∂r = v X ∂ ∂r X + v Y ∂ ∂r Y + v Z ∂ ∂r Z , r r 3 ∂ ∂v = r X r 3 ∂ ∂v X + r Y r 3 ∂ ∂v Y + r Z r 3 ∂ ∂v Z s ∂ ∂v = s X ∂ ∂v X + s Y ∂ ∂v Y + s Z ∂ ∂v Z
F sr , F sv , F rv , F vr , F vv -F rr .

dim Lie(F 0 , F 1 , ...) = 5 < dim R 6 = 6

Moreover, rank deficiency implies that an integral of motion exists, which happens to be the

L F 1 (det(r, v, ŝ)) = s X ∂ ∂v X det(r, v, ŝ) + s Y ∂ ∂v Y det(r, v, ŝ) + s Z ∂ ∂v Z det(r, v, ŝ) = s X (-r Y s Z + r Z s Y ) + s Y (r X s Z -r Z s X ) + s Z (-r X s Y + r Y s X ) = 0

  Control sets for different reflectivity coefficients and s = 1. Here, uX is the projection of u towards ŝ, while uY and uZ are orthogonal components.

Fig. 1 :

 1 Fig. 1: Components of the SRP force.

Fig. 2 :

 2 Fig. 2: Euler angles orienting the orbit with respect to the reference frame. Here, h and e denote

Fig. 3 :

 3 Fig. 3: Schematic representation of a half-space of the neighborhood of I where motion is (locally)

Fig. 4 :

 4 Fig. 4: Approximation of the control set (blue) by a convex cone (red).

Fig. 5 :

 5 Fig. 5: Example of the solution of Problems (8) (black curve) and (9) (red dot). Here,

  i.e., a, b H = Re(a), Re(b) + Im(a), Im(b) , and denote Φ

9 Fig. 8 :

 98 Fig. 8: Results for different geocentric orbits

  γ3 fixed at 120[deg]

Fig. 9 :

 9 Fig. 9: Minimum cone angle as a function of Euler angles.

F- 1 r 3 FF

 3 denoting ŝ • r a scalar product of two vectors ŝ and r, computation of Lie brackets gives the following results: rv ;F 1 = s ∂ ∂v = F sv ; F 01 = F 0 , F 1 = -F sr F 001 = F 0 , [F 0 , F 1 ] = 3 (ŝ • r) r 5 F rv -F sv r 3 , F 101 = F 1 , [F 0 , F 1 ] = 0 0001 = F 0 , [F 0 , [F 0 , F 1 ]] = 1 r 3 F sr + 3(v • r) r 5 F sv + 3(ŝ • v) r 5 -15(ŝ • r)(v • r) r 7 F rv + 3(ŝ • r) r 5 (F vv -F rr )

  projection of the angular momentum h towards ŝ, namely ŝ• h = ŝ • (r × v) = det(r, v, ŝ) In fact,Lie derivative of ŝ • h with respect to the controlled vector field,F 1 = s X ∂ ∂v X + s Y ∂ ∂v Y + s Z ∂ ∂v Z is:

Clearly, because of the periodicity of the free motion, the angular position determined by the true anomaly does

In practice, of course, the check can only be made at a single point, while the condition must hold for all I to ensure global controllability.

We note that Gkl = G(-k,--l) because G(I, f )u is real valued.

of square functional systems for trigonometric polynomials was proposed to efficiently evaluate this requirement. Extensive exploitation of the algorithm revealed that any orbital condition can tolerate a certain amount of absorption of solar radiation. A remarkable byproduct of this analysis is that local controllability properties hold universally (that is for all orbital conditions) for non-ideal sails provided a sufficient amount of incoming radiation is reflected.

VII. Acknowledgements

This work was partially supported by ESA (contract no. 4000134950/21/NL/GLC/my).

APPENDIX A: LIE BRACKETS COMPUTATION FOR A PERFECTLY-ABSORPTIVE

SAIL

Consider a control-affine dynamical system

where M is an n-dimensional manifold, F i : M → T M are smooth vector fields on M. We note that Eq. ( 2) can be recast into the form of Eq. (A1) by choosing columns of G(I, f )R(I, f ) as vector fields F i .

To compute Lie algebra, only directions of vector fields matter. To simplify, assume that control system is given by a simple two-body equation with a term of perturbation, as denoted in (A2) with ŝ solar vector, considered fixed for a few orbits. For a perfectly absorptive solar sail, only the cross-sectional surface is controlled, so that control is assumed to be u ∈ [0, 1] with certain coefficient ε defining SRP magnitude.

Using x = (r, v) ∈ R 6 as state vector, where r and v denote Cartesian position and velocity