María Teresa Becerra
email: mtbecerra@gmail.com

Jesús Mavárez
email: jesus.mavarez@univ-grenoble-alpes.fr

BECERRA and MAVAREZ: A NEW ESPELETIA FROM COLOMBIA Espeletia saboyana (Millerieae, Asteraceae), a New Critically Endangered Caulirosula from Cordillera Oriental, Colombia

m tall, with whitish-cinereous appearance, sessile leaves with elliptic laminae, cymose capitulescences with aphyllous unbranched portion and 3-5 capitula, 5-13 phyllaries, the external ones 17.4-26.0 mm long, 49-92 ray florets 13.0-18.4 mm long, and 162-347 disc florets 8.2-12.2 mm long. Similar to E. incana Cuatrec., but with longer laminae, bigger foliar area, and fewer ray floret series. Also similar to E. murilloi Cuatrec., but with longer ray tubes, ray paleae, disc corollas, disc styles, disc anther appendages, and disc paleae. In addition, E. saboyana is distributed allopatrically in regards with E. incana and E. murilloi. This species is proposed to be listed as IUCN Critically Endangered (CR), since it is restricted to a single population with Extent of Occurrence = 1.83 km 2 and Area of Occupancy = 0.48 km 2 . Besides, this population is severely fragmented into 40-50 demes, the largest one covering only 34,000 m 2 . Most demes are found in humid soils and wetlands submitted to continuous habitat deterioration as water is drained for use in agriculture. Espeletia saboyana is rare or absent in these newly drained areas, which are either colonized by other species that thrive on drier soils or transformed into pastures for cattle grazing. Furthermore, analyses of size class distributions show that recruitment rates of E. saboyana are significantly lower in small demes, which, if not reversed, can further reduce population size in the near future.

Keywords-Andes, Biological Conservation, Caulescent rosette, Demography, High-elevation, Morphometrics, Páramo, Taxonomy.

The genus Espeletia Mutis ex Bonpl. is endemic to the high elevations in the Andes of Venezuela, Colombia and Northern Ecuador, with some outliers in Venezuelan Cordillera de la Costa and Colombian Sierra Nevada de Santa Marta.

These plants are well known for their extraordinary diversity of growth-forms, which include profusely branched, dichotomous and unbranched trees, shrubs, branched caulirosulas, sessile caulirosulas and, particularly, the emblematic giant caulirosulas. Espeletia currently includes 139 species, usually found from ca. 2400 m above sea level (hereafter masl) in the upper Andean forest to ca. 4700 masl near the historical glacier limits, although most species occur between 3200-4000 masl, i.e. within the neotropical high-elevation open-vegetation ecosystem known as páramo [START_REF] Diazgranados | A nomenclator for the frailejones (Espeletiinae Cuatrec., Asteraceae)[END_REF][START_REF] Cuatrecasas | A Systematic Study of the Subtribe Espeletiinae (Heliantheae, Asteraceae)[END_REF][START_REF] Mavárez | A taxonomic revision of Espeletia (Asteraceae). The Venezuelan radiation[END_REF], 2021). Espeletia species frequently become ecologically dominant elements of páramos, usually reaching adult densities > 1 plant/m 2 (Smith 1981, Fagua and[START_REF] Fagua | Growth rates, reproductive phenology and pollination ecology of Espeletia grandiflora (Asteraceae), a giant andean caulescent rosette[END_REF]. Thus, many species, and notably the giant caulirosulas, contribute so much to the characteristic physiognomy of páramos in terms of abundance and biomass that for the general public these plants represent the defining element of this beloved ecosystem. This fact is celebrated in numerous postcards, books, webpages, and even Colombian and Venezuelan national banknotes.

Taxonomic research should be nearly complete in Espeletia after more than 200 years since its formal description [START_REF] Humboldt | Guidelines for Using the IUCN Red List Categories and Criteria[END_REF], and especially when considering the morphological distinctiveness of the genus, its ecological dominance, and the great interest it generates to the people. But in reality much taxonomic research remains to be done in Espeletia. For instance, several species are known by only a handful of specimens, while many localities are severely under-represented in botanical collections, especially in Colombia, where access to some páramos has been difficult due to the lack of roads and decades of armed conflict. In fact, the taxonomic knowledge of Colombian Espeletia is on constant change as numerous species are continuously being described from this country. For instance, 24 out of the 26 Espeletia species described since 2000 are from Colombia (Díaz-Piedrahita and Pedraza 2001, Díaz-Piedrahita and Obando 2004, Díaz-Piedrahita et al. 2006, Díaz-Piedrahita and Rodríguez-Cabeza 2008[START_REF] Zar | Biostatistical Analysis 5 th edition[END_REF][START_REF] Diazgranados | A new species of Coespeletia (Asteraceae, Millerieae) from Venezuela[END_REF][START_REF] Diazgranados | A new species of Espeletiopsis (Millerieae, Asteraceae) from Colombia[END_REF][START_REF] Diazgranados | Espeletia praesidentis, a new species of Espeletiinae (Millerieae, Asteraceae) from northeastern Colombia[END_REF][START_REF] Mavárez | Taxonomic novelties in páramo plants. Espeletia ramosa (Asteraceae), a new species from Colombia[END_REF][START_REF] Alzate | Espeletia restricta (Millerieae, Asteraceae), a new species from the páramos of northern Colombia[END_REF][START_REF] Becerra | Espeletia ocetana (Millerieae, Asteraceae), a New Tall Caulirosula from Colombia[END_REF]. This trend probably reflects a delay in taxonomic research that is being rapidly remedied with improvements in road infrastructure and better personal safety conditions in Colombia.

As part of an ongoing research program on the diversity and taxonomy of

Espeletia, we recently had the opportunity to sample a population of giant caulirosulas with sessile leaves and aphyllous capitulescences from a páramo west of the town Saboyá, Boyacá department, Colombia (Fig. 1). This population, hereafter referred to as Saboyá Espeletia, has been known for quite some time but its taxonomic determination has remained rather enigmatic. For instance, [START_REF] Rangel | Estrategia Corporativa para la Caracterización con Fines de Manejo y Conservación de Áreas de Páramo en el Territorio CAR[END_REF] and [START_REF] Arellano | Patrones en la distribución de la vegetación en áreas de páramo de Colombia: heterogeneidad y dependencia espacial[END_REF] could not assign it to any known species and entered it as Espeletia sp. in their analyses. On the other hand, [START_REF] Rojas | Plan de Conservación y Manejo de los Frailejones del Territorio CAR[END_REF] determined that Saboyá Espeletia belongs to E. incana Cuatrec., which is a somewhat morphologically similar species that nonetheless seems to differ in indumentum color, leaf size, leaf venation structure, and capitulescence size. Besides, the nearest confirmed E. incana population is located in Páramo de la Rusia, about 80 km E-NE of Saboyá. An evaluation of the other giant caulirosulas with aphyllous capitulescences from Western Boyacá suggests that Saboyá Espeletia is also morphologically similar to E. murilloi Cuatrec., a species with known populations in Páramo de Merchán, about 15 km E of Saboyá. However, as above, E. murilloi also appears to show some differences in regards with the size and shape of its leaves and capitulescences.

To test whether Saboyá Espeletia is conspecific with E. incana or E. murilloi we performed uni-and multivariate morphometric analyses of several leaf, capitulescence, and capitulum traits with great taxonomic importance in Espeletia (Cuatrecasas 2013, Becerra and[START_REF] Mavárez | A taxonomic revision of Espeletia (Asteraceae). II. Updated list of taxa, nomenclature, and conservation status in the Colombian radiation[END_REF]. Our results clearly show that the Saboyá population does not belong to either of the two mentioned species. The Espeletia from Páramo de Saboyá is described and illustrated in this work. Morphological traits-From each individual we sampled two leaves on opposite sides of the middle part of the rosette (without sheaths to avoid excessive damage) and two fully developed capitulescences. Most leaf and capitulescence traits were measured directly on fresh samples without magnification using a tape or a ruler. Leaf lamina area was estimated indirectly, as the ratio between the weight of a 1:1 leaf drawing measured on a 0.01 g precision balance and the paper density in g/cm 2 . In addition, we sampled and preserved in 70% ethanol one mature capitulum from each individual of the three taxa. Capitulum traits were measured in the laboratory using a 3.5 magnification device and a caliper. The ray floret series was estimated as the number of complete whorls plus the fraction of incomplete ones (e.g. 3.25, 4.50). Three replicates per capitulum were obtained for traits associated with phyllaries, florets, and paleae. In total, 43 quantitative traits were measured on each plant: nine in leaves, six in capitulescences, and 28 in capitula, with replicates per individual varying between one (11 capitulum traits), two (all leaf and capitulescence traits), and three (17 capitulum traits) (Table 1).

MATERIALS AND METHODS

Study sites and sampling-We

Morphological analyses-Statistical analyses followed the procedures described in [START_REF] Becerra | Espeletia ocetana (Millerieae, Asteraceae), a New Tall Caulirosula from Colombia[END_REF]. In short, the data for each trait were averaged per individual, log10 transformed, and their normality and homocedasticity verified using Shapiro-Wilk and Levene tests, respectively (no deviations found). Univariate morphological differences among taxa were evaluated using one-way analysis of variance (ANOVA), with significance cut-off values adjusted for the number of tests using a Bonferroni procedure ( adj = 0.05/43 = 1.16e-3, Sokal & Rohlf 2012, p239), whereas post-hoc comparisons for significant tests were analyzed using Dunnett's tests considering Saboyá Espeletia as control (Zar 2010, p234). Multivariate morphological differences among taxa were evaluated with Principal Component Analysis (PCA). Both uni-and multivariate tests and analyses were performed in R using the basic package Stats and visualized with ggplot2 v3.3.2 and factoextra v1.0.7.

Taxonomic descriptions-Complete leaves (with sheaths), mature capitulescences, and capitula of ten additional individuals of Saboyá Espeletia were sampled and described following the terminology employed in [START_REF] Roque | in Systematics, Evolution, and Biogeography of Compositae[END_REF] and [START_REF] Becerra | Espeletia ocetana (Millerieae, Asteraceae), a New Tall Caulirosula from Colombia[END_REF]. Most descriptions were based on dry material, although some stem and rosette traits were measured directly in the field, while descriptions and drawings of capitulum traits were made from samples preserved in 70% ethanol. These specimens were deposited in herbarium COL, whereas duplicates will be sent to FMB and UPTC (herbaria acronyms follow Thiers 2021).

Distribution and conservation-To know the geographic distribution of

Saboyá Espeletia, we inspected satellite images in and around Páramo de Saboyá using Google Earth, searching for areas above 3000 masl with open vegetation that could represent suitable habitats for Espeletia. Favorable areas were visited in June 2021 and any Espeletia species found was collected and identified. We followed IUCN guidelines to estimate Extent of Occurrence (EOO), i.e. the area of the smallest polygon containing all presence sites and in which no internal angle exceeds 180 degrees, a.k.a. minimum convex polygon (IUCN 2019), and Area of Occupancy (AOO), i.e. the sum of the area of all 0.1 km wide square cells in which the taxon is present. EOO and AOO were calculated with GeoCat [START_REF] Bachman | Supporting Red List threat assessments with GeoCAT: geospatial conservation assessment tool[END_REF], using the geographic coordinates recorded by GPS of all patches with confirmed individuals of Saboyá Espeletia.

To obtain some preliminary insights on the demographic functioning and short-term future trends of Saboyá Espeletia, we conducted a field study describing the spatial variation in population density and size-classes of this taxon. 26 plots (10-15 m × 4 m) were established in twelve patches or demes chosen according to their size: eleven plots were placed within four "large" patches (8218-33,644 m 2 ) and fifteen plots within eight "small" patches (724-3364 m 2 ) (Fig. 2). We counted all individuals in each plot, measured their stem height, and determined their reproductive status. The smallest flowering plant (i.e. smallest adult) had a stem of 0.95 m, so 1 m was designated as the minimum size at maturity. To facilitate data analysis each individual was assigned to one of the following size classes defined according to reproductive status and stem height: juvenile-1: < 0.5 m, juvenile-2: 0.5-1.0 m, adult-1: 1.0-1.5 m, and adult-2: > 1.5 m. The differences in abundance of size classes between large and small patches were evaluated using a nested ANOVA, with factor "plot" (levels: 1 to 26) nested within factor "patch size" (levels: small or large).

RESULTS

ANOVA-Results of F-tests from one-way ANOVA are shown in Table 2. 76%

of the 86 pairwise comparisons between Saboyá Espeletia and either of the other two species showed statistically significant differences. The magnitudes of the observed differences are very high, as half of the 65 significant pairwise comparisons showed adjusted p-values < 1e-10 and only two had values > 1e-4. In regards to trait type, 78% and 79% of pairwise comparisons involving leaf and capitulum traits showed significant differences, respectively, whereas only 58% of differences in capitulescence traits were significant. Indeed, only five out of 43 traits analyzed (i.e. LWL, CTPL, HPN, HRN, and HDN) showed no statistically significant differences between Saboyá Espeletia and neither of the compared species. Saboyá Espeletia differed significantly from E. incana and E. murilloi at seven and seven out of nine leaf traits analyzed, respectively, with no-overlapping values for lamina area (LLA) and length (LLL), and nearly so for max width (LLW), midpoint width (LMW), and base width (LBW). Differences in capitulescences were somewhat less important, as Saboyá Espeletia differed from E. incana and E. murilloi at four and three out of the six traits analyzed, respectively, all of them showing overlapping values among the three taxa. Conversely, Saboyá Espeletia differed significantly from E. incana and E. murilloi at 20 and 24 of the 28 capitulum traits analyzed, respectively, with non-overlapping values for the number of ray floret series (HRS), and the lengths of ray floret tubes (HRTL), ray and disc paleae (HRPa1L, HRPa2L, HDPaL), disc florets (HDL), disc styles (HDSL), disc stamens (HDStL), and disc anthers (HDAL).

PCA-A graphical representation of the PCA results describing the ensemble of differences between and within the three taxa is shown in Fig. 3A-B. The first two axes in this analysis, explaining 61.7% of morphological variance among taxa, clearly show that Saboyá Espeletia individuals segregate into a distinct morphologic cluster without overlapping with E. incana or E. murilloi clusters. The average within cluster Euclidean distance in PCA space for Saboyá individuals is 1.59 (sd. 1.04), which is 4.7 and 6.4 times smaller than the average distances against centroids of E. incana (7.41, sd. 1.57) and E. murilloi (10.23, sd. 1.69), respectively. Furthermore, without exception, each Saboyá individual is closer to its own centroid than to those of the other species, and its nearest neighbor is always another member of the same population. On the other hand, the morphological distinctiveness of Saboyá Espeletia against E. incana or E. murilloi also holds when PCA are built using leaf and capitulescence traits (Fig. 3C) or capitulum traits (Fig. 3D).

In sum, the ensemble of ANOVA and PCA results reveal an important morphological discontinuity between Saboyá Espeletia and either E. incana or E. Diagnosis-Caulescent rosette up to 5 m tall, with whitish-cinereous indumentum; sessile leaves with elliptic laminae; cymose capitulescences with aphyllous unbranched portion and 3-5 capitula; 5-13 phyllaries, the external ones 17.4-26.0 mm long; 49-92 ray florets, 13.0-18.4 mm long; 162-347 disc florets, 8.2-12.2 mm long. Similar to E. incana, but with longer laminae (31.0-48.5 cm vs.

19.0-30.5 cm), larger foliar area (147-389 cm 2 vs. 59-130 cm 2 ), and fewer ray floret series (2.6-3.9 vs. 4.1-6.3). Also similar to E. murilloi, but with longer ray tubes (3.2-7.0 mm vs. 1.7-2.9 mm), longer outmost ray paleae (10.8-16.2 mm vs. 6.8-10.7 mm), longer disc corollas (8.2-12.0 mm vs. 6.2-8.1 mm), longer disc styles (8.5-10.9 mm vs. 6.2-8.3 mm), longer disc anther appendages (2.6-4.8 mm vs. 1.8-2.5 mm), and longer disc paleae (8.2-11.9 mm vs. 6.0-8.1 mm).

Caulirosula (Fig. 4A-B), polycarpic, unbranched stem, erect, up to 5 m tall, covered with marcescent leaves, diam about 8-10 cm below the rosette base.

Rosette with whitish-cinereous indumentum, diam 80-90 cm, height 50-70 cm.

Leaves sessile (Fig. 4C,5A); sheath open, trapezoid, plurinerved, moderately coriaceous, 6.5-8.9 cm × 6.0-8.3 cm, abaxially floccose with golden trichomes 20-25 mm long, adaxially with golden trichomes 12-15 mm long only near the lamina base; lamina elliptic, coriaceous but flexible, 31.0-49.0 cm × 6.0-12.5 cm, length/width ratio 3.2-5.7:1, gradually attenuate towards the base, width at base 1.2-2.8 cm, apex acute, margins entire, undulate; abaxial surface with dense indumentum lanate-floccose, green-yellowish, trichomes < 15 mm long near the lamina base, < 5 mm long elsewhere, midvein pale-yellow, round and prominent all along, secondary veins prominent, parallel, mostly craspedromous, distally curved or divided, 8-14 veins/dm at lamina mid-section, angles 56°-80°; adaxial indumentum dense, lanate-floccose, soft, whitish-cinereous, trichomes < 25 mm long, shorter and sparser in oldest leaves, midvein visible all along, slightly prominent and gradually broadened towards the lamina base, secondary veins hardly visible. Capitulescences cymose (Fig. 4E-G, 5B-C), simple, opposite, axillary, erect, 6-12 coetaneous, axes robust, rigid, 24-59 cm long, capitulescence/lamina ratio 0.6-1.5, covered with abundant floccose whitish-yellowish indumentum, trichomes 10-15 mm long; proximal (unbranched) section 19-48 cm long, aphyllous; distal (branched) section 5-16 cm long, with 3-5 capitula, usually as 3 in a terminal cyme, and the others, if present, on 2 monocephalous lower branchlets, occasionally as 1 terminal capitulum and 2 monocephalous lower branchlets; central pedicels 1.0-4.5 cm long, lateral pedicels of terminal cymes 2.4-4.0 cm, sometimes sessile, pedicels of lower branchlets 3.5-7.0 long; subtending bracts opposite, sheaths connate and amplectant, lanceolate, acuminate, both sides covered with whitish lanate-floccose indumentum, abaxial trichomes < 10 mm long, adaxial trichomes < 15 mm long, bracts on lower branchlets 7.0-9.0 cm × 1.5-2.2 cm, bracts on 3-capitulum cymes 3.0-7.5 cm × 1.4-3.0 cm, bracts on single terminal capitula 2.4-4.7 cm × 1.1-2.1 cm. Capitula radiate (Fig. 4D), usually nodding, each with 222-426 flowers, involucre diam 17.4-31.3 mm, ray circle diam 25-44 mm, disc diam 16-34 mm; involucre campanulate, phyllaries 5-12(13) (Fig. 6A,6E), in two rows, outer phyllaries (2)3-4, herbaceous, 16.2-30.8 mm × 9.4-19.4 mm, ovate-oblong, apex acute, abaxially densely barbate, adaxially glabrous proximally, distally barbate, inner phyllaries 2-8(10), herbaceous, 12.0-21.1 mm × 5.5-13.5 mm, oblong, apex subobtuse, abaxially villous, adaxially villous only at apex; first series of ray floret paleae 13-22 (Fig. 6B), oblong, thin, herbaceous, subamplectant, 10.8-16.2 mm × 3.7-8.0 mm, subacute, abaxially villous along central longitudinal line and apex, adaxially glabrous; second and third series of ray floret paleae oblong (Fig. 6C), membranaceous, whitish-hyaline, amplectant, 9.8-13.4 mm × 2.1-4.3 mm, 3-5veined, central vein thicker, abaxially barbate at apex, trichomes erect, yellowish, < 1.2 mm long, adaxially glabrous; disc floret paleae oblong (Fig. 6D), membranaceous, whitish-hyaline, amplectant, 8.2-11.9 mm × 1.2-2.9 mm, 3-5veined, abaxially barbate distally, trichomes erect, yellowish-brown, < 1.2 mm long, adaxially glabrous; receptacle slightly convex, diam 10.5-18.9 mm. Ray florets 49-92 (Fig. 6F), in 2.6-3.9 series, 12.6-18.7 mm long; corollas liguliform, yellow, 8.9-14.9 mm long; tubes 3.2-7.0 mm long, with scattered trichomes < 0.5 mm long, eglandular, erect, transparent, uniseriate, septate, with round apex; limbs oblong, 4.9-8.8 mm 1.4-2.7 mm, 2-3-dentate, 5-7-veined, with scattered sessile glands; styles 4.9-10.9 mm long, branches 1.7-2.7 mm long; cypselae 3-ribbed (Fig, 6G), 2.8-5.0 mm × 1.5-2.8 mm, ellipsoidal-prismatic, dark brown, epappose, three denticulate at apex. Disc florets 162-347 (Fig. 6H), 8.2-11.1 mm long; tubes 2.4-4.7 mm, with scattered trichomes in the distal end, < 0.5 mm long, eglandular, uniseriate, septate, with round apex; limbs tubular-infundibuliform, 4.6-6.6 mm long, lobes triangular < 0.7 mm long, margins thickened, glabrous; stamens 7.2-11.5 mm long, anthers brown, 2.6-4.8 mm long, exceeding the corolla, calcarate, base shortly sagittate, apical appendages ovate-acute; styles 7.5-12.1 mm long, arms short, conical, glabrous, dorsally papillose, nectaries tubular, 5-denticulate, 1.0-1.5 mm long. juvenile-2 size class was found in the two plots analyzed in large patch L4 (Fig. 7B), which is located on a seemingly former humid area that was drained and is now densely colonized by grasses. This unusual patch was excluded in further analyses.

Etymology-Espeletia

The remaining data showed important differences in regards with the density of size classes between large and small patches (Fig. 7A). For instance, size class distribution in large patches exhibits the well-known "inverted J" shape, in which juvenile individuals account for approximately 66% of the population, while the juvenile-1 class alone represents nearly 47%. This trend is completely reversed in small patches, where a "J" shaped size class distribution was found, with adults accounting for 66% of the population and adult-2 being by far the most abundant class (52%). ANOVA analyses show that differences in size class abundances between large and small patches of E. saboyana are highly significant (juvenile-1: p = 1.94e-3, juvenile-2: p = 0.083, adult-1: p = 0.019, adult-2: p = 1.16e-5).

IUCN Red List Category-Critically

Endangered (CR) (criteria B1ab(i-iii, v)+B2ab(i-iii, v).
Espeletia saboyana is restricted to a single population in the Colombian Andes with Extent of Occurrence EOO = 1.83 km 2 and Area of Occupancy AOO = 0.48 km 2 , estimated from 347 presence records with GeoCat. However, the area actually occupied by E. saboyana, calculated directly from the sum of the areas of its inhabited polygons is only 0.165 km 2 . In fact, the only population of this species is severely fragmented into 40-50 patches or demes, the largest of which covering only 34,000 m 2 (Fig. 7B). The majority of these patches are located on humid soils and are surrounded either by areas of relatively natural habitat with drier soils in which E. saboyana is very rare (i.e. demes < 10 adults inds., Fig. 8B), and that are instead dominated by the rosette Espeletia argentea, dense stands of the invasive alien species Ulex europaeus (Fig. 8A,8E), and some grasses, or by areas impacted by agricultural activities in which E. saboyana is completely absent (Fig. 8D). There is continuous deterioration in the quality of the habitat of E. saboyana, since water from humid areas in Páramo de Saboyá is drained for use in agriculture (Fig. 8E). Newly drained areas are either colonized by other species that thrive on drier soils or by pastures for cattle grazing. Furthermore, demographic analyses of size class distributions show that recruitment rates of E. saboyana are significantly lower in small patches, which, if not reversed, can further reduce their size, and eventually make them disappear. incana, respectively (Table 2), with some of the overlaps observed at other traits due to a single extreme value in either species. The PCA results also unquestionably demonstrate that E. saboyana represents a distinct and non-overlapping morphological cluster compared to E. incana and E. murilloi, which provides strong support to the recognition of E. saboyana as a separate species. Furthermore, PCA results show that the average morphological distance between E. saboyana and the other species is several times larger than the average distance among E. saboyana individuals, revealing very important morphological discontinuities between taxa.

Additional Specimens Examined-Espeletia

All these results also hold for PC analyses performed with vegetative or reproductive morphological traits, separately. Thus, this study provides another convincing example of the potential of morphometric analyses to distinguish closely related taxa in Asteraceae (e.g. [START_REF] Owen | A multivariate morphometric analysis of the Symphyotrichum boreale -S. nahanniense -S. welshii complex (Asteraceae: Astereae)[END_REF][START_REF] Cook | A multivariate morphometric analysis of Solidago subsect. Glomeruliflorae (Asteraceae: Astereae)[END_REF][START_REF] Semple | A multivariate morphometric analysis of the Solidago uliginosa complex (Asteraceae: Astereae; S. subsect. Maritimae) in eastern North America[END_REF], and the second one in Espeletia [START_REF] Becerra | Espeletia ocetana (Millerieae, Asteraceae), a New Tall Caulirosula from Colombia[END_REF].

Espeletia saboyana belongs to the group of tall caulirosulas with aphyllous and opposite capitulescences, which is one of the most diverse in Espeletia, particularly in Colombia [START_REF] Cuatrecasas | A Systematic Study of the Subtribe Espeletiinae (Heliantheae, Asteraceae)[END_REF]. However, E. saboyana can be easily distinguished from other morphologically similar species distributed around Páramo de Saboyá (i.e. P. Merchán, P. Rabanal, P. Iguaque), using simple combinations of vegetative and reproductive traits. For instance, it differs from E. incana by having longer laminae (> 31 cm vs. < 31 cm), larger foliar area (> 147 cm 2 vs. < 130 cm 2 ), and fewer ray floret series (< 3.9 vs. > 4.1). It also differs from another somehow similar species, E. murilloi, by the bigger lengths of its ray tubes (> 3.2 mm vs. < 2.9 mm), outmost ray paleae (> 10.8 mm vs. < 10.7 mm), disc corollas (> 8.2 mm vs. < 8.1 mm), disc styles (> 8.5 vs. < 8.3 mm), disc anther appendages (> 2.6 mm vs. < 2.5 mm), and disc paleae (> 8.2 mm vs. < 8.1 mm). The only other tall caulirosula with aphyllous capitulescenses present around Páramo de Saboyá is E. barclayana Cuatrec., a species that can easily be distinguished from E. saboyana by its pseudopetiolate leaves and its sub-corymboid capitulescences with [START_REF] Cuatrecasas | A Systematic Study of the Subtribe Espeletiinae (Heliantheae, Asteraceae)[END_REF][START_REF] Mavárez | A taxonomic revision of Espeletia (Asteraceae). II. Updated list of taxa, nomenclature, and conservation status in the Colombian radiation[END_REF]).

The addition of E. saboyana rises to 140 the number of recognized species in the genus Espeletia, with 86 in the Colombian or southern clade [START_REF] Mavárez | A taxonomic revision of Espeletia (Asteraceae). II. Updated list of taxa, nomenclature, and conservation status in the Colombian radiation[END_REF]) and 54 in the Venezuelan or northern clade [START_REF] Mavárez | A taxonomic revision of Espeletia (Asteraceae). The Venezuelan radiation[END_REF]. According to current records there are 92 Espeletia species distributed within Colombia's political borders.

Conservation-Espeletia saboyana is endemic to the highest elevations in a single, small, and geographically isolated páramo in the Colombian Cordillera Oriental. Its Extent of Occurrence EOO = 1.83 km 2 and Area of Occupancy AOO = 0.48 km 2 are both below the maximum values for the IUCN CR criteria, i.e. B1: EOO < 100 km 2 and B2: AOO < 10 km 2 . Moreover, there is (a) only one population of this species, composed itself by a series of rather isolated patches (Fig. 2), and (b) continuing projected decline in (i) EOO, (ii) AOO, and (iii) quality of habitat; as the wetlands that constitute the main habitat of E. saboyana are being drained by cutting ditches into the ground to collect and transport water for agriculture (Fig. 8E). In fact, we noticed that this species is very rare in drained areas, which are instead inhabited by another Espeletia species with preference for dry soils (i.e. E. argentea), by an aggressive alien species that can form large mono-specific stands in degraded páramos (i.e. Ulex europaeus) (Barrera-Cataño et al. 2019, Figs. 8A, 8E), and by grasses used as feeding grounds for cattle (Fig. 8D). Furthermore, juveniles are very rare or even absent in isolated, small or dried patches of E. saboyana (Fig. 7), probably due to some negative fitness consequences associated with reduced habitat size (e.g. [START_REF] Fahrig | Effects of Habitat Fragmentation on Biodiversity[END_REF] or small deme size (e.g. [START_REF] Aguilar | Plant reproductive susceptibility to habitat fragmentation: review and synthesis through a metaanalysis[END_REF][START_REF] Krushelnycky | Evaluating the interacting influences of pollination, seed predation, invasive species and isolation on reproductive success in a threatened alpine plant[END_REF]). However, the absence of juveniles in large patch L4, an area that was formerly humid and is now dry and entirely colonized by grasses, suggests that E. saboyana can already suffer the negative consequences of reduced habitat quality before the subsequent reduction in patch size. In any case, regardless the causes, current demographic functioning in E. saboyana allows inferring (v) a projected continuous decline in the number of mature individuals of this species in the near future. In summary, the information available makes it possible to affirm undoubtedly that this species should be registered as Critically Endangered CR, according to the IUCN criteria B1ab(i-iii, v)+B2ab(i-iii, v) (IUCN 2019). Espeletia saboyana joins therefore the unfortunate list of critically endangered Espeletia, which in Colombia now includes 17 out of 84 species with proficient data, i.e. excluding seven Data Deficient DD species [START_REF] Becerra | Espeletia ocetana (Millerieae, Asteraceae), a New Tall Caulirosula from Colombia[END_REF]Mavárez 2021, iucnredlist.org).

The discovery of this new species and its categorization as critically endangered illustrate the importance of detailed taxonomic knowledge while evaluating the conservation status of a species, particularly for taxa with restricted distributions. The Colombian regional environmental agency with jurisdiction over Páramo de Saboyá did not include E. saboyana in its list of Espeletia species with priority for conservation since it was considered as a marginal population of the more broadly distributed E. incana [START_REF] Rojas | Plan de Conservación y Manejo de los Frailejones del Territorio CAR[END_REF]. Furthermore, in this particular case the negative consequences of taxonomic misidentification may work as well in the opposite direction, by inflating the distribution area of E. incana. Indeed, some tall caulirosulas with whitish indumentum in Boyacá dept., e.g. E. pescana (S.Díaz) S.Díaz, are frequently misidentified as belonging to E. incana. The real distribution and number of localities of E. incana are probably much smaller than currently assessed in IUCN analyses (Vulnerable category, VU), and this species would very likely need to be reclassified as Endangered (EN) [START_REF] Mavárez | A taxonomic revision of Espeletia (Asteraceae). II. Updated list of taxa, nomenclature, and conservation status in the Colombian radiation[END_REF].

According to the evidence above, threats shed upon E. saboyana do not concern the species itself (e.g. trade, consumption, direct use for human activities), but instead affect its habitat. The survival of this species is therefore highly dependent on the conservation of the extent and quality of the páramo habitat in Páramo de Saboyá. In line with Colombian national regulations that recognize páramos and wetlands as areas of special ecological importance (Presidencia de la República de Colombia 2010, Congreso de la República de Colombia 2018), we recommend the following actions: (i) to study in more detail the population dynamics of E. saboyana, with special attention on the causes of reduced rates of recruitment/juvenile survival in small patches, (ii) to study the potential impact of introgressive hybridization with E. argentea upon E. saboyana's genotypes and phenotypes, (iii) to evaluate the conservation status of wetlands in Páramo de Saboyá, and accordingly formulate and implement measures to rapidly restore the natural extent and quality of theses habitats, (iv) to implement awareness raising activities involving local communities and other relevant stakeholder in pro of the conservation of the species and its habitat, and (iv) to support the implementation of effective monitoring systems aimed at evaluate the outcomes and success of these and other conservation actions [START_REF] Robinson | How to ensure threatened species monitoring leads to threatened species conservation[END_REF]. 3.0/6.0 (1.6)/9.8 (<5e-10) 7.0/10.9 (2.1)/15.0 7.8/12.5 (3.8)/24.0 (n.s.) CTPL 1.3/2.6 (1.1)/5.6 (n.s.) 1.5/3.9 (1.9)/7.6 1.5/6.2 (3.3)/12.6 (n.s.) CHN 3.0/3.1 (0.2)/4.0 (<6e-10) 3.0/4.8 (0.5)/5.0 3.0/4.8 (1.5)/8.5 (n.s.) HRD 2.4/2.8 (0.2)/3.1 (<4e-8) 2.6/3.4 (0.4)/4.2 2.4/2.8 (0.3)/3.4 (<1e-6) HID 1.7/1.9 (0.2)/2.1 (<1e-7) 1.8/2.4 (0.3)/3.1 1.6/2.0 (0.3)/2.4 (<5e-5) HDD 1.5/1.7 (0.1)/1.9 (<3e-9) 1.6/2.2 (0.3)/2.7 1.4/1.8 (0.2)/2.2 (<2e-7) HReD 11.5/14.3 (1.8)/18.6 (n.s.) 10.5/14.2 (2.2)/18.9 8.0/11.3 (1.7)/14.1 (<1e-5) HReH 3.9/5.5 (0.8)/7.2 (n.s.) 4.6/6.4 (1.1)/8.1 3.2/4.2 (0.7)/5.9 (<1e-10) HPN 6/8.8 (2.1)/15 (n.s.) 5/7.9 (2.3)/13 5/7.7 (2.0)/12 (n.

HRS

4.1/4.9 (0.6)/6.3 (<1e-10) 2.6/3.6 (0.3)/3.9 2.9/4.1 (0.5)/4.8 (<9e-4) HRL 12.4/14.3 (1.2)/16.4 (n.s.) 13.0/15.1 (1.5)/18.4 8.4/10.4 (1.4)/13.8 (<1e-10) HRLL 5.9/7.4 (1.0)/9.6 (n.s.) 5.2/6.5 (0.9)/8.7 3.2/5.1 (1.2)/8.0 (<1e-5) HRLTN 2.0/2.3 (0.3)/2.7 (<3e-5) 2.3/2.8 (0.3)/3.3 2.0/2.6 (0.4)/3.3 (n.s.) HRTL 2.6/4.0 (0.5)/4.5 (<2e-7) 3.6/4.8 (0.8)/6.6 1.8/2.3 (0.3)/2.7 (<1e-10) HRCL 3.0/3.3 (0.2)/4.0 (<1e-7) 3.2/3.9 (0.4)/4.9 2.7/3.1 (0.2)/3.5 (<1e-10) HRSL 4.5/6.4 (0.6)/7.1 (<3e-3) 5.4/7.2 (0.8)/8.5 4.0/4.8 (0.5)/5.9 (<1e-10) HRBL 2.2/2.6 (0.2)/3.0 (<1e-5) 1.9/2.2 (0.2)/2.6 1.2/1.7 (0.2)/1.9 (<1e-10) HRPa1N 6/9.9 (3.0)/15 (<1e-10) 13/17.1 (2.9)/22 10/13.7 (2.9)/20 (<4e-3) HRPa1L 9.3/11.0 (0.9)/12.1 (<2e-8) 11.5/13.0 (0.8)/14.6 7.2/8.5 (0.8)/9.9 (<1e-10) HRPa2L 7.5/9.3 (0.9)/11.3 (<1e-10) 10.3/11.4 (0.8)/13.0 6.9/7.6 (0.6)/8.9 (<1e-10) HDN 157/228.5 (60.1)/331 (n.s.) 162/255.0 (52.2)/347 109/214.9 (52.2)/297 (n.s.) HDL 7.6/8.7 (0.7)/9.9 (<1e-4) 8.5/9.5 (0.7)/10.8 6.5/7.1 (0.4)/7.9 (<1e-10)

  collected 20 flowering Saboyá Espeletia individuals on August 2019 and September 2020 at 3350 masl in Páramo de Saboyá, also known as Páramo de Telecom (5.714881°N, 73.815718°W), municipality of Saboyá, Boyacá department, Colombia (Fig. 1, filled square). 20 flowering individuals of E. incana were collected on November 2020 and November 2021 at 3650 masl in Páramo de la Rusia (5.929764°N, 73.083890°W), in the western side of the Guantiva-la Rusia páramo complex, Boyacá, Colombia (Fig. 1, filled circle). 20 flowering individuals of E. murilloi were collected as follows (Fig. 1, filled triangles): three on October 2020 at 3486 masl in Páramo La Cortadera (5.551906032°N, 73.10549296°W), five on October 2020 at 3150 masl in Páramo Siscunsí (5.53229497°N, 72.72864603°W), six on June 2021 at 3383 masl in Páramo de Merchán (5.704074958°N, 73.65707399°W), five on July 2021 at 2926 masl near Arcabuco (5.729302°N, 73.387125°W), and one on November 2020 at 3435 masl in Páramo Rabanal (5.438630981°N, 73.54459397°W).

  murilloi, indicating that they are not co-specific. Therefore, the individuals from Saboyá are regarded as a different Espeletia species whose description is provided below. TAXONOMIC TREATMENT Espeletia saboyana M.T. Becerra & Mavárez, sp. nov. TYPE: COLOMBIA. Boyacá dept., Saboyá, Páramo de Saboyá, in the road Saboyá-Jesús María, 3350 masl, 5.714881°N, 73.815718°W, August 19 2019, M.T. Becerra, J. Mavárez & J. Aguirre 53 (holotype: COL, isotypes: COL, FMB, UPTC).

  saboyana specific epithet refers to the Páramo and Municipality of Saboyá, the only locality where the species is known to occur. Distribution and Habitat-Espeletia saboyana is only distributed in the central and highest elevation section of Páramo de Saboyá, Municipio de Saboyá, Departamento de Boyacá, Colombia. It is absent from the northern and southern sections of Páramo de Saboyá, as well as from nearby Páramo de Merchán. This species grows between 3220-3350 masl, on flat areas and gentle slopes within páramo and sub-páramo habitats with humid to very humid soils rich in organic matter, in which it coexists with Espeletia argentea Bonpl., Blechnum loxense (Kunth) Hieron., Chusquea tessellata Munro, Bucquetia glutinosa DC., Ulex europaeus L., Hypericum sp. and Puya sp, and that are surrounded by secondary-growth forests with Pernettya prostrata Sleumer, Linochilus rosmarinifolium Benth., Ageratina sp. and Baccharis sp. Demography-The mean population density of E. saboyana, based on data derived from 26 plots covering 1248 m 2 , is 0.97 (sd. 0.22) individuals/m 2 . The four size classes represent 32.4% (juvenile-1), 14.2% (juvenile-2), 17.3% (adult-1), and 36.2% (adult-2) of the 1221 sampled individuals. No individual of juvenile-1 or

  saboyana-Colombia, BOYACÁ: Saboyá, road Saboyá-Jesús María, 3350 masl, 5.714881°N, 73.815718°W, M.T. Becerra, J. Mavárez & J. Aguirre 54 (Paratype: COL, FMB, UPTC); M.T. Becerra, J. Mavárez & J. Aguirre 55 (Paratype: COL, FMB, UPTC); Páramo de Saboyá, 3000 masl, Stancik 1061 (COL); Páramo Cerro de Telecom, 3250 masl, Camargo 7269 (COL); Páramo de Telecom, Veredas Molino y Tibista, [5.712778°N, 73.816944°W], 3370 masl, Santos-C. 75 (COL). -Espeletia murilloi-Colombia, BOYACÁ: Páramo de Arcabuco, 2800 masl, E. Pérez-Arbeláez & J. Cuatrecasas 8098 (Isotype: COL); Páramo de la Rusia, R. Romero-Castañeda & R. Romero-Mejía 6812 (COL); S. Díaz-Piedrahita 38 (COL); Santuario de Fauna y Flora de Iguaque, J.P. Vargas & J. Hernández JPVN-20 (COL); Duitama, 3300-3500 masl, B.V. Rodríguez-Cabeza BVR-1935 (COL). -Espeletia incana-Colombia, SANTANDER: Páramo de la Rusia, western slope, 3300-3500 masl, J. Cuatrecasas 10430 (Isotype: COL); BOYACÁ: Duitama, Páramo de la Rusia, 3570 masl, A. M. Cleef 7150 (COL); Hoya de la Quebrada de Becerra, 3500-3600 masl, J. Cuatrecasas & R. Jaramillo-Mejía 28728 (COL). DISCUSSION Analyses of uni-and multivariate datasets comprising 43 leaf, capitulescence, and capitulum traits clearly support the recognition of E. saboyana as a morphologically distinct species easily distinguishable from E. incana and E. murilloi. The ANOVA results demonstrate that E. saboyana is significantly different from at least one of the other two species at 39 out of 44 traits analyzed, with the majority of significant comparisons showing extremely low p-values. Nine and four traits showed no range value overlap between E. saboyana and E. murilloi or E.

  disc florets (HDN) (n) 1 Number of capitula (CHN) (n) 3 Disc floret length (HDL) (mm) 1 Capitulum 3 Disc floret limb length (HDLL) (mm) 1 Involucre diameter (HID) (cm) 3 Disc floret tube length (HDTL) (mm) 1 Ray circle diameter (HRD) (cm) 3 Disc floret style length (HDSL) (mm) 1 Disc diameter (HDD) (cm) 3 Disc floret stamen length (HDStL) (mm) 1 Receptacle diameter (HReD) (mm) 3 Disc floret anther appendage length (HDAL) (mm) 1 Receptacle height (HReH) (mm) 3 Disc palea length (HDPaL) (mm)

FIG. 1 .

 1 FIG. 1. Maps of Colombian Cordillera Oriental showing known localities of Saboyá

Figure 1

 1 Figure 1

Figure 2

 2 Figure 2

Figure 3

 3 Figure 3

Figure 4

 4 Figure 4

Figure 5

 5 Figure 5

Figure 6

 6 Figure 6

Figure 7

 7 Figure 7

Figure 8

 8 Figure 8

  > 5 relatively large capitula. All other Espeletia species in this region have opposite and thyrsoid capitulescences with foliose unbranched portions, i.e. E. argentea, E. boyacensis Cuatrec., E. cayetana (Cuatrec.) Cuatrec., E. chocontana Cuatrec., E. grandiflora, E. raquirensis Rodr.-Cabeza & S.Díaz, E. tunjana Cuatrec., or alternate and corymboid capitulescences, i.e. E. corymbosa Bonpl., E. garciae Cuatrec., E. glandulosa Cuatrec., E. guacharaca S.Díaz, E. muiska Cuatrec., E. pleiochasia Cuatrec., E. rabanalensis S.Díaz & Rodr.-Cabeza

TABLE 1 .

 1 List of traits analyzed and units in the three Espeletia populations 591 evaluated. Superscripts in plant parts and traits represent replicates per individual 592 and per trait, respectively. 593

	2 Leaf lamina	1 Number of phyllaries (HPN) (n)
	1 Area (LLA) (cm 2 )	3 External phyllaries in row 1 length (HEP1L) (mm)
	1 Length (LLL) (cm)	3 External phyllaries in row 2 length (HEP2L) (mm)
	1 Max. width (LLW) (cm)	1 Number of ray florets (HRN) (n)
	1 Length/width ratio (LLR)	1 Number of ray floret series (HRS) (n)
	1 Length at max. width (LWL) (cm)	3 Ray floret length (HRL) (mm)
	1 Midpoint width (LMW) (cm)	3 Ray floret limb length (HRLL) (mm)
	1 Lamina base width (LBW) (mm)	1 Number of ray floret limb teeth (HRLTN) (n)
	1 Veins at midsection (LVD) (n/dm)	3 Ray floret tube length (HRTL) (mm)
	1 Vein angles at midsection (LVA) (°) 3 Ray floret cypsela length (HRCL) (mm)
	2 Capitulescence	3 Ray floret style length (HRSL) (mm)
	1 Length (CL) (cm)	3 Ray floret style branch length (HRBL) (mm)
	1 Proximal section length (CPL) (cm) 1 Number of ray paleae in 1st row (HRPa1N) (n)
	1 Distal section length (CDL) (cm)	3 Ray palea in 1st row length (HRPa1L) (mm)
	1 Capitulescence/lamina ratio (CLR) 3 Ray palea in 2nd row length (HRPa2L) (mm)
	1 Terminal pedicels length (CTPL)

TABLE 2 .

 2 Minimum, mean (sd), and maximum values of morphological traits averaged per individual in E. incana, Saboyá Espeletia and E. murilloi. Superscripts indicate P-values from Dunnett's tests against Saboyá Espeletia. Grey background indicates no-overlap between Saboyá Espeletia and the other species. Trait acronyms given in Table 1.

			Min/Mean (sd.)/Max
		E. incana	Saboyá Espeletia	E. murilloi
	LLA	65.9/90.1 (18.6)/129.8 (<1e-10) 151.9/222.8 (49.2)/382.5 51.0/98.2 (24.2)/149.1 (<1e-10)
	LLL	21.0/26.6 (3.1)/30.5 (<1e-10)	31.3/37.1 (3.1)/44.8	23.3/28.1 (2.9)/33.5 (<1e-10)
	LLW	4.4/6.4 (1.0)/8.1 (<1e-8)	6.2/9.0 (1.3)/12.4	3.6/5.0 (0.9)/6.9 (<1e-10)
	LLR	3.2/4.3 (0.8)/6.8 (n.s.)	3.5/4.2 (0.5)/5.3	4.2/5.8 (1.1)/8.3 (<2e-7)
	LWL	49.0/53.9 (4.6)/64.0 (n.s.)	52.3/57.9 (3.7)/64.5	47.4/54.7 (3.9)/65.0 (n.s.)
	LMW	4.2/6.3 (1.1)/8.1 (<4e-7)	6.0/8.5 (1.3)/12.2	3.4/4.8 (0.9)/6.5 (<1e-10)
	LBW	1.1/1.4 (0.3)/2.0 (<1e-4)	1.3/1.8 (0.3)/2.5	0.6/1.0 (0.2)/1.4 (<1e-10)
	LVD	11.0/15.0 (2.3)/19 (<9e-8)	8.0/11.2 (1.6)/14	10.3/14.5 (1.9)/17 (<1e-6)
	LVA	42.0/54.4 (5.3)/63.0 (<1e-9)	59.0/66.6 (4.0)/75.5	54.7/63.9 (6.0)/73.8 (n.s.)
	CL	22.0/28.3 (4.3)/40.2 (<1e-10)	30.0/43.1 (6.0)/52.7	41.5/54.5 (9.7)/77.3 (<3e-5)
	CLR	0.8/1.1 (0.2)/1.4 (n.s.)	0.8/1.2 (0.1)/1.5	1.4/1.9 (0.3)/2.6 (<1e-10)
	CPL	14.5/22.3 (5.1)/36.0 (<8e-7)	23.0/32.2 (5.6)/43.8	26.3/42.0 (9.5)/67.3 (<5e-4)
	CDL

ACKNOWLEDGMENTS

We are very grateful for the help provided by Julián Aguirre from Universidad Nacional de Colombia, Catherine Numa from IUCN, John Alexander Sánchez and Olga Lucía Hernández from Instituto Humboldt, and María Isabel Castro from Universidad de La Salle Colombia. We also thank Nicolás Contreras, Juana Muñoz, Diego Ramírez, Edilberto Suarez and Carlos Suárez for help during fieldwork, and the staff at Instituto Alexander von Humboldt and herbarium COL for help with plant preservation and lab facilities. Collections were granted by the Colombian Autoridad Nacional de Licencias Ambientales (ANLA), Permiso Marco de Recolección de Especímenes, Resolución 01435 (September 3/2018, attributed to J. Aguirre) and Resolución 01194 (July 13/2020, attributed to M. T. Becerra). This work is part of the Ph.D. dissertation of M. T. Becerra.

AUTHOR CONTRIBUTIONS

MTB provided morphological data, demographic data, line drawings, photographs, and was the primary author of the new species description. JM provided maps and performed statistical analyses. Both authors contributed to the text.


domain: spi phys

text:
Lucas Fuster
email: lucas.fuster@laplace.univ-tlse.fr

Gerjan Hagelaar

Romain Pascaud

Antoine Simon

Patrick Hoffmann

Laurent Liard

Olivier Pascal

Thierry Callegari

Plasma-based microwave power limitation in a printed transmission line: a self-consistent model compared with experimental data

Plasma-based microwave power limitation in a suspended microstrip transmission line integrating a micro hollow cathode discharge (MHCD) in its center is experimentally and numerically studied. Transient and steady state microwave power measurements exhibit a limitation threshold of 28 dBm and time responses of 25 microseconds. Intensified charge-coupled device (ICCD) imaging shows that microwave breakdown occurs at the top of the MHCD. The plasma then extends towards the microwave source within the suspended microstrip transmission line. Besides, a self-consistent model is proposed to simulate the non-linear interaction between microwave and plasma. It gives numerical results in agreement with the measurements, and show that the plasma expansion during the transient response is related to a shift between the ionization source term and the electron density maximum. The propagation speed, under the tested conditions, depends mainly on the stepwise ionization from the excited states.

Introduction

With the development of high-power microwave (HPM) sources [1], the susceptibility of microwave receivers becomes an important issue [2]. As a result, microwave power limiters are sometimes introduced in receivers before sensitive components to avoid permanent damage [3]. A good microwave power limiter should have low insertion loss below a given input power threshold so as not to increase the noise figure of the receiver, and above this limiting threshold it should quickly reduce its output power, known as leakage power, below a safe value.

Microelectromechanical systems (MEMS) [4], high-temperature superconductors [5,6], or phasetransition materials [7] have been used as active elements for microwave power limitation, but most current solutions use solid-state components such as positive-intrinsic-negative (PIN) diodes [8][9][10]. PIN diodes generally have modest insertion loss and fast response times while offering small dimensions, relatively simple implementation, and low cost. However, their power handling may not be sufficient with regard to the HPM sources currently under development [11].

Because of their high power handling capability, plasma-based microwave power limiters are also investigated [START_REF] Smullin | Microwave Duplexers[END_REF][START_REF] Patel | Microstrip plasma limiter IEEE MTT[END_REF][14][15][16][17][START_REF] Simon | Static and dynamic control of limiting threshold in plasma-based microstrip microwave power limiter[END_REF][19][20][21][22]. For such limiters, the incident HPM signal causes gas breakdown, and this microwavedriven plasma discharge then absorbs and/or reflects the incoming power, thus limiting the output power. Plasma-based limiters were first implemented in waveguide technology in transmit-receive (TR) tubes, which are still used in duplexers for high-power pulsed radar systems [START_REF] Smullin | Microwave Duplexers[END_REF]. However, with dimensions of the order of several centimeters, TR tubes remain bulky and not suitable for integration on printed circuit boards, the most widespread technology for the design of microwave receivers.

As far as we know, Patel et al. designed the first plasma-based microwave power limiter in planar printed circuit technology, but its high limiting threshold of 47 dBm clearly restricted its use [START_REF] Patel | Microstrip plasma limiter IEEE MTT[END_REF]. Since then, this limiting threshold has been successfully reduced and/or controlled using resonant circuits [14], pre-ionization techniques [20], or both [15,[START_REF] Simon | Static and dynamic control of limiting threshold in plasma-based microstrip microwave power limiter[END_REF]. In resonant circuits, the magnitude of the electric field is enhanced at resonance, thus reducing the microwave power required for gas breakdown [14,[START_REF] Simon | Static and dynamic control of limiting threshold in plasma-based microstrip microwave power limiter[END_REF]. Regarding pre-ionization, it provides seed electrons such that the breakdown electric field drastically decreases [20,23]. It is interesting to note that similar solutions have also been implemented in two-dimensional (2D) photonic crystals (PhC) which are investigated for the realization of millimeter-wave receivers [17,19,21].

Although these plasma-based limiters have interesting performances, most of current work is based on experimental studies and few on numerical models that may provide interesting information on the interaction between microwave power and plasma discharge. For instance, numerical diagnostics of the electron density have been proposed [16,17,21]. In that case, the plasma is modelled as a uniform three-dimensional (3D) object described by the Drude model in an electromagnetic software, and its electron density is sought in simulation in order to obtain the same response as in measurement. Recently, Gregório et al. have proposed to couple fluid and Maxwell equations to study the interaction between microwave and plasma [24]. This model was applied to the simulation of a PhC power limiter, but no comparison with measurements was provided. Finally, the lack of validated models that faithfully reproduce the non-linear interaction within plasma-based microwave power limiters makes their design and optimization difficult.

In this paper, we present the design of a plasmabased microwave power limiter in suspended microstrip technology. In addition to having good performance as a limiter, its geometry allows it to be studied with a self-consistent 2D numerical model.

A direct comparison between experimental and numerical results shows that this model is predictive and that it can be used for studying the proposed limiter.

Device under test and experimental setup

Device under test

The exploded view of the proposed plasma-based microwave power limiter is shown in figure 1b. It is similar to the one studied by Simon et al. in [20], but the microstrip transmission line has been replaced by a suspended microstrip transmission line as described below.

The microwave power limiter in suspended microstrip technology, also called device under test (DUT) hereafter, is made up of several planar layers stacked on top of each other. As shown in figure 1b, the suspended microstrip line consists of two conductors, namely a ground plane and a line printed on a thin dielectric substrate. Dielectric spacers are placed on both ends of the transmission line to provide a constant gap between the two conductors, and thus form a low pressure gas-filled cavity that hosts the plasma discharge during the microwave power limitation. Note that the cavity remains open on its sides to allow optical measurements of the plasma during the experiments. Due to the dielectric spacers, the input and output of the whole transmission line are similar to microstrip transmission lines, whereas the middle part is a suspended one. The dielectric substrates are Rogers RO4003C with ε r = 3.55 and tan δ = 0.0021. The thickness of the dielectric spacers is h gap = 1.524 mm, their length is L sup = 50 mm, and the thickness of the suspended dielectric substrate is h sms = 0.508 mm. The width of the printed line varies according to the dielectric configuration between the conductors in order to achieve a 50 Ω characteristic impedance along the whole transmission line, and thus allows the impedance matching of the DUT. Therefore, we have W line1 = 4.5 mm for the microstrip transmission line section and W line2 = 8 mm for the suspended one. Finally, the microwave part of the DUT is W sms = 100 mm wide, L sms = 200 mm long, and h tot = h gap + h sms = 2.032 mm thick, whereas the cavity is 100 mm long. Note that the transmission line is connected at both ends to 50 Ω coaxial straight connectors.

The principle of operation of this plasma-based microwave power limiter has already been exposed in [15,20], and it is here briefly recalled. When a microwave signal is propagating along the transmission line, namely along x direction, the electric field is mostly confined between the ground plane and suspended microstrip and oriented along y direction. If the incoming microwave power is large enough, for instance during an HPM attack, gas breakdown may be ignited in the low pressure gas-filled cavity. This microwave-driven plasma then absorbs and/or reflects the HPM signal thereby reducing the output power. In order to reduce the microwave power limiting threshold, pre-ionization is used [20,21,23]. It consists in generating seed electrons inside the cavity with a DC micro-hollow cathode discharge (MHCD) [25].

The DC circuit for pre-ionization, that is to say the MHCD, uses a cylindrical copper cathode printed on a Rogers RO4003C dielectric substrate below the ground plane. Its diameter is Φ cathode = 80 mm. As observed in figure 1b, the ground plane of the transmission line is used as the anode for the MHCD. A cylindrical hole is drilled in the center of the cathode through the dielectric substrate and the two electrodes. The dimensions of the MHCD hole have been chosen to h M HCD = 1.524 mm and Φ M HCD = 1 mm to lower the breakdown voltage at 10 Torr in Argon. The MHCD is thus in normal glow regime with DC current I DC = 3 mA [26].

The DUT was measured with a vector network analyzer (Rohde & Schwarz, ZVL13) and simulated using Ansys HFSS, a 3D frequency-domain electromagnetic solver. Figure 2 presents the magnitude of the smallsignal reflection and transmission coefficients (i.e., S 11 and S 21 parameters, respectively) of the DUT without pre-ionization. A good agreement is observed with a measured reflection coefficient lower than -20 dB and transmission coefficient better than -0.4 dB for the bandwidth of the experimental setup, namely 2 to 4 GHz.

Experimental setup

Figure 3 shows the experimental setup. The DUT is placed in a quartz vacuum chamber evacuated to a base pressure of 10 -5 Torr by a vacuum pump station (Edwards Vacuum, TSM1D1001) and then backfilled with Argon to the desired pressure, namely 10 Torr, which is measured by a piezo transducer gauge (Pfeiffer Vacuum, CMR361). The Argon flow is fixed to 0.1 L/min using a mass flow controller (Bronkhorst, EL-FLOW Prestige). The DUT is thus totally immersed in a large volume of Argon, including the cavity between the ground plane and the suspended microstrip line. A high voltage generator (Technix, SR2KV-2KW) is connected to the MHCD electrodes through a DC vacuum feedthrough. It is used to ignite and sustain the DC plasma discharge for preionization. The current I DC flowing through the discharge is limited by a series resistor R = 10 kΩ. Regarding the high-power microwave measurement devices, they allow for both steady-state and transient measurements.

The steady-state measurement setup involves the blue and black components in figure 3. A continuous wave (CW) microwave signal is generated by a microwave generator (Anapico, APSIN20G) and then amplified by a power amplifier (MC2 Technologies, HPA60W2-4G). This power amplifier (PA) has a maximum output power of 60 W between 2 and 4 GHz. The HPM signal is then delivered to the DUT through 50 Ω coaxial vacuum feedthroughs. An isolator (Quest Microwave, NM2040C02) has been inserted between the PA and the first coupler in order to protect the microwave amplifier against high reflected power. Finally, leakage power of the DUT is dissipated by a high-power 50 Ω termination (Weinschel Associates, WA1448). As shown in figure 3, the input, reflected, and output powers of the DUT are coupled by means of two -30 dB directional couplers (SHX, DDTO-3-30).

One is connected between the PA and the input port of the DUT while the other is connected between its output and the 50 Ω termination. Power measurements are done thanks to three power sensors (Keysight, U8481A) connected to the coupled ports of the two couplers. Insertion loss from coaxial cables, feedthroughs, and couplers have been measured with a vector network analyzer (Rohde & Schwarz ZVL13). All steady-state power measurements reported thereafter have thus been corrected so that they refer to DUT ports. Note that the frequency range is limited to 2 -4 GHz due to the PA, isolator, and couplers characteristics, and the maximum power at the DUT input is 45 dBm because of coaxial feedthroughs power handling.

For transient measurements, slight modifications are done to the experimental setup (red and black components in figure 3).

A microwave switch (American Microwave Corporation, SWMDJV-1DT-2ATT) is first inserted between the microwave signal generator and the PA. The resulting signal is a highpower microwave pulse. In addition to that, the power sensors are replaced by an oscilloscope (Keysight, MSO9254A) that records the input, reflected, and output signals of the DUT. Finally, a camera (Princeton Instruments, PI-MAX 512) with a red-blue optimized intensifier (Princeton Instruments, GEN II) is added to the experimental setup in side view of the DUT. The microwave switch, the oscilloscope, and the camera are synchronized with a control signal whose pulse width and repetition frequency are equal to 1 ms and 100 Hz, respectively. Each picture results from a 200 ms exposure time with a 20 ns gate time.

Model and numerical setup

In order to study the non-linear response of the plasmabased microwave power limiter described in figure 1, a 2D self-consistent numerical model was implemented. This model has already been detailed in [27], but we remind here its principles and point out the specific adaptations we made for the present study.

As the plasma density and chemistry evolve on a long time scale (microsecond) compared to that of the microwave oscillations (sub-nanosecond), our model is composed of two sub-models describing these two time scales separately. Both sub-models are based on a combination of fluid equations and Maxwell equations. The sub-model for the plasma evolution is obtained by averaging these equations over the short (microwave) time scale, while that for the microwaves neglects the time variations of the plasma parameters. The submodels are coupled via the absorbed microwave power (activating the plasma), and the electron density and collision frequency (affecting the microwaves). Thus, during the simulation of the plasma evolution, the microwave model is periodically called to update the microwave power absorption profile, using the latest plasma profile.

In the following sections, we first describe the model geometry, then the two sub-models, and finally the plasma-chemical processes that are taken into account.

Simulated geometry

The suspended microstrip transmission line shown in figure 1 has a large width-to-height ratio w line2 /(h gap + h sms ) of 3.9. As a result, fringing effects can be neglected and the central longitudinal part of the DUT can be seen as an inhomogeneously dielectric filled parallel-plate transmission line [START_REF] Collin | Foundations for Microwave Engineering[END_REF]. We thus propose to simulate it using the 2D simulation domain described in figure 4 that represents a part in the (x, y) plane of the suspended microstrip transmission line. This 2D model assumes that the geometry is invariant in the z direction, and that plasma and microwave characteristics are uniform along z, but still functions of x and y. Since the dominant mode is the TM x mode, E x , E y , and B z field components are considered. Note that this mode behaves as a quasi-TEM mode in the low frequency limit, without plasma, and that no higher order modes are excited as long as k 0 (h gap + h sms ) ≪ 1, with k 0 the wave number in vacuum [START_REF] Collin | Foundations for Microwave Engineering[END_REF]. This condition is satisfied at 3 GHz, our frequency of interest.

In order to reduce the computation time, the length of the plasma cavity region to the right of the MHCD source (x > 70 mm in figure 4) has been shortened to 9 mm. We verified that it does not modify the numerical results for our microwave power range. The dielectric constant of the two dielectric spacers at the edges of the device has been set equal to 1 in the simulations, instead of 3.55 in the experiment. This modification allows us to reproduce the impedance matching (when there is no plasma in the cavity) that is experimentally obtained by changing the width of the printed line of the DUT.

Finally, as shown in figure 4, the microwave field is excited by an artificial current sheet near the left edge of the simulation domain as explained in section 3.2. Regarding the MHCD, it is not simulated in detail but modelled by adding a localized source term in plasma equations which reproduces a realistic electron density and ensures that the gas breakdown occurs at the same incident power as in experiments. This MHCD source term has a Gaussian profile centered at x = 70 mm in figure 4, near the ground plane, with a width of 1 mm in the x direction, and a width of 0.3 mm in the y direction. The amplitude of the source term is of 1.35 × 10 23 m -3 .s -1 .

Microwave model

The microwaves are described by solving the Maxwell equations [27]:

∇ × Ẽ Ẽ Ẽ = - ∂ B B B ∂t (1) 1 µ 0 ∇ × B B B = J J Jexc -en e ṽ ṽ ṽe + ε 0 ε r ∂ Ẽ Ẽ Ẽ ∂t (2)
where Ẽ Ẽ Ẽ and B B B are the oscillating electric and magnetic fields, n e the electron density obtained from the plasma model, ṽ ṽ ṽe the oscillating component of the electron mean velocity, J J Jexc an external oscillating current which excites the microwaves and µ 0 is the magnetic permeability of vacuum. The second term on the right hand side of equation ( 2) is an approximation of the plasma current density, neglecting the ion contribution (because m i ≫ m e ) and oscillations of the electron density (ñ e ≪ n e ).

The oscillating electron velocity is obtained from the local momentum equation ∂ ṽ ṽ ṽe ∂t = -e m e Ẽ Ẽ Ẽν e ṽ ṽ ṽe ,

where ν e is the momentum transfer frequency of collisions with the gas. This set of equations is solved by the finite-difference time-domain (FDTD) method, using a Yee mesh that coincides with the mesh of the plasma model, with perfect conductor boundary conditions at the metal boundaries (top and bottom) and first-order Mur absorption boundary conditions at the open boundaries (left and right) [START_REF] Taflove | Computational Electrodynamics: The Finite-Difference Time-Domain Method[END_REF].

For the excitation current J J Jexc in equation ( 2) we use a thin current sheet near the left edge of the simulation domain (see figure 4). This generates an electric field in the y direction which then propagates both to the right and to the left, while inducing other field components and interfering with the wave reflected by the plasma. We fix the amplitude of J J Jexc at an arbitrary value, hence the amplitudes of Ẽ Ẽ Ẽ, B B B, and ṽ ṽ ṽe are also arbitrary and need to be normalized afterwards. This is done by imposing a fixed incident power, as follows.

First the amplitude of the reflected wave is determined by combining the net fields in two planes P 1 and P 2 located at a small distance δ to the left and to the right of the current sheet, respectively (see figure 4)

E 2 ref l = 1 2 sin 2 (βδ) ( Ẽy,2 -Ẽy,1 ) 2 t ( 4
)
B 2 ref l = 1 2 cos 2 (βδ) ( Bz,1 + Bz,2 ) 2 t ( 5
)
where the angle brackets ⟨.⟩ t indicate the time-average over the microwave period and β is the propagation phase constant of the dominant mode within this transmission line, given in the low frequency limit by [START_REF] Collin | Foundations for Microwave Engineering[END_REF] β

≈ k 0 h gap + h sms h gap + h sms /ε r (6)
that is to say β ≈ 1.1k 0 for our problem.

The time-averaged reflected power is then found by integrating the reflected Poynting vector over the cross section of the transmission line

P ref l = αL z 2µ 0 E ref l B ref l dy (7)
where L z is the effective size of the model geometry in the z direction (≈ width of the transmission line) and α is a normalization factor to correct for the arbitrary amplitude of the excitation current. Finally, subtracting (7) from the net power passing through the P 2 , we get an expression for the incident power which allows us to determine the normalization factor

α 1 α = L z µ 0 P in ⟨ Ẽy,2 Bz,2 ⟩ t - 1 2 E ref l B ref l dy, ( 8
)
where P in is a constant incident power, one of the main input parameters of our simulations.

Once the normalization factor is known, the microwave power absorbed by the plasma follows from P abs = L z P e (x, y)n e (x, y)dxdy (9) where the integrals are over the plasma region and

P e = α -e Ẽ E E • ṽ v v e t ( 10
)
is the time-averaged absorbed power per electron, which is injected into the energy equation (??) of the plasma model. In addition, the transmitted power P out is calculated by integrating the Poynting vector over the plane P 3 at the right domain edge. Our simulations accurately reproduce steady-state energy conservation P in = P ref l + P abs + P out .

For optimal consistency, we calibrated the value of L z using full 3D simulations of the DUT without plasma with Ansys HFSS. From these simulations we determined the magnitude of the electric field Ẽ3D at the position of the MHCD as a function of the input power P 3D , and then adjusted L z such as to reproduce these results

µ 0 ⟨ Ẽ2 ⟩ t L z ⟨ Ẽy Bz ⟩ t dy = ⟨ Ẽ2 3D ⟩ t P 3D . (11)
From this calibration, we found a value of L z = 7.34 mm, which we then fixed in the simulations. (The geometrical value is 8 mm, similar but a little different.)

Chemical model

The electrons e and five different Argon species are considered in our simulations: atomic ions Ar + , molecular ions Ar + 2 , two excited species Ar * m and Ar * r , and background gas Ar. The complete set of chemical reactions between these species is presented in Table 1.

The rate coefficients of reactions 1 to 9 were calculated from cross section data with the BOLSIG+ software [41], which solves the electron Boltzmann equation in order to take into account the shape of the electron energy distribution function (EEDF) selfconsistently. These rate coefficients are used in the form of lookup tables as a function of T e , assumed to be equal to 2/3 times the electron mean energy.

The excited Argon states are grouped into two species: an effective metastable species Ar * m and an effective resonant species Ar * r . These species are assumed to consist mainly of the first four excited Argon states, namely the 1s 5 and 1s 3 metastable states and the 1s 4 and 1s 2 resonant states (in Paschen notation), respectively. However, the rate coefficients for the production of these species (reactions 2 and 3) take into account not only direct excitation of the 1s states but also cascading from higher states, using the effective excitation cross sections of the Morgan database [START_REF]Morgan database[END_REF].

The rate coefficients for stepwise ionization (reactions 6 and 7) are based on cross sections for each of the 1s states, calculated using the analytical solution proposed by Deutsch et al. [START_REF] Deutsch | [END_REF]32] which is consistent with experimental measurements in Argon [33].

The rate coefficients of reactions 8 and 9 were calculated from cross sections for electron-impact state-to-state transitions between the different 1s states taken from the detailed BSR database [START_REF]BSR database (Quantum-mechanical calculations[END_REF].

The radiative decay rate of A * r (reaction 10) was estimated as the average of the natural decay rates of the 1s 4 and 1s 2 states, approximately 3.32 × 10 8 s -1 [START_REF]NIST database[END_REF], multiplied by an escape factor to account for the effect of resonance radiation imprisonment. For an infinite slab geometry and line broadening due to collisions, the effective escape factor is (see equations (5.17) in [START_REF] Holstein | [END_REF] and (6.5) in [36])

g ef f = 1.150 π λ ef f 3h gap (12)
where λ ef f = 105.8 nm is the average wavelength of the resonance radiation.

Results and discussion

As mentioned previously, both simulations and measurements of the DUT were carried out considering Argon at 10 Torr. In addition, a frequency of 3 GHz was chosen for the HPM signal. All the following results were obtained for these conditions. Unless otherwise stated, the MHCD is always switched on with a DC current I DC equal to 3 mA.

Transient results

Plasma dynamics and power ratios evolution

The DUT was first measured using the experimental setup of figure 3 in transient configuration. Optical measurements were achieved with the fast camera to observe the dynamics of the plasma. Both simulation and measurements were conducted using a 40 dBm HPM pulse. The time t = 0 µs subsequently refers to the start of the HPM pulse. Note that unlike our experimental setup, the simulation gives us access to the electron density of the plasma as a function of time.

In figure 5, we can see for instance the measured light emitted by the plasma as well as the simulated electron density taken at different times. When microwave gas breakdown occurs in the device, the plasma grows and then moves towards the left, from the position of the MHCD toward the input of the DUT. For an input power P in = 40 dBm, the maximum value of the simulated electron density is about 10 19 m -3 .

In addition to that, we measured the instantaneous reflected and output powers, denoted as P ref l and P out , respectively. Experimentally, the absorbed power is computed as P abs = P in -P out -P ref l . Formally, this is the definition of the dissipated power, but here, the radiated power from the transmission line remains low and the dissipated power can reasonably be considered as mainly absorbed by the plasma discharge. The simulated reflected, absorbed and output powers have also been measured. From these power results, we can compute the reflected, output, and absorbed power ratios relative to the input power that are denoted T ref l , T out , and T abs , respectively. Figure 6 shows the time evolution of the different power ratios during microwave breakdown with a good agreement between simulation and measurements, even if a discrepancy of 10%, that remains to be understood, appears around 100 µs. We observe the typical response of a mainly absorbing microwave power limiter with its output power ratio that is decreasing over time up to its leakage value while its dissipated power ratio is highly increasing. The power of the HPM signal at the output of the DUT is here reduced by 90 % after 25 µs. Beyond 150 µs, the power ratios no longer change.

The time evolution of the power ratios may be linked to the plasma dynamics. Therefore, in addition to the images of figure 5, figure 7 shows the time evolution of the average electron density over the plasma domain ⟨n e ⟩ x,y as well as the maximum of the electron density max(n e ). One can distinguish four steps in the plasma dynamics (the times given are those of the simulation):

(1) Before the beginning of the HPM pulse, the MHCD is in steady-state and it produces a localized electron density of up to 9 × 10 17 m -3 . From t = 0 µs to t 1 = 2.31 µs, there is the microwave breakdown with an increasing electron density that reaches a maximum at t 1 (see figure 7). This corresponds to the strong increase of the absorbed power ratio that reaches a local maximum at t 1 . As shown in figure 5, this maximum of electron density at this stage remains at the location of the MHCD, that is to say of the pre-ionization region until t 2 = 10.9 µs, and the absorbed power ratio remains approximately constant. (2) From t 2 to t 3 = 16.11 µs, the maximum of electron density propagates on the right hand side of the MHCD and starts to decrease. Simultaneously, we observe an increase of the electron density on the left hand side. At t 3 , the two peaks of electron density have the same amplitude as shown in figure 5. (3) From t 3 to t 4 = 26.70 µs the electron density finishes to decrease on the right hand side and to increase on the left hand side of the MHCD location. The absorbed power ratio thus increases and the reflection and transmission coefficients both decrease. At t 4 the mean electron density over the plasma domain ⟨n e ⟩ x,y reaches its constant value, as shown in figure 7. (4) Finally, as shown in figure 5, this pattern propagates from the location of the MHCD towards the input of the DUT and a steady-state is reached.

Propagation mechanisms

In this section, we propose to give a more precise description of the physical mechanisms involved in the propagation of the pattern. The numerical code allows us to look at the relative importance of each reaction in the production of electrons. With the reaction rates obtained with the EEDF self-calculated by BOLSIG+ in the local field approximation, no electrons are created by the direct ionization while 61.1 % and 33.2 % are produced by the ionisation of metastable and resonant states, respectively. These values have been obtained by integration over the plasma domain.

In figure 8, we present the electron density profile, the excited states density profile and the source terms of the electrons and the excited states due to the excitation and the stepwise ionisation, defined as:

S + Ar * m,r = (k 2 + k 3 )n Ar n e (13)
S + e = (k 6 n Ar * m + k 7 n Ar * r )n e (14) where k 2 , k 3 , k 6 and k 7 are the rates of reaction 2, 3, 6 and 7 (see table 1). All the profiles have been averaged over the plasma domain along the y direction. They have been computed at P in = 40 dBm and t = 75.2 µs. The excited states appear in the first half of the electron density profile, because of the absorption of the electromagnetic energy by the electrons. The source term for the electrons ⟨S + e ⟩ y reaches a maximum at x = 46.8 mm, whereas the electron density profile ⟨n e ⟩ y reaches a maximum at x = 47.8 mm. It is this gap between the position of the maxima of ⟨n e ⟩ y and ⟨S + e ⟩ y that induces the pattern movement. Note that the plasma does not propagate by transport of ions and electrons, but rather by an ionization wave that propagates along the cavity. For input powers between 29 dBm and 42 dBm, the front velocities have been extracted both numerically and experimentally and are presented in figure 9.

The mechanism of the pattern propagation may also be described analytically. According to [42], a fixed point stability analysis for the following nonlinear diffusion equation:

∂n ∂t -D ∂ 2 n ∂x 2 = ν i n -r ei n 2 (15)
with D the diffusion coefficient, ν i the ionisation frequency, and r ei the recombination coefficient, all of them being uniform, allows to conclude that the plasma front propagates with a velocity of the order of v ∼ 2 √ Dν i . Boeuf et al. have shown that the diffusion coefficient is of the form [43]:

D ef f ≈ α α + 1 D e + 1 α + 1 D a (16)
with D e and D a the free and ambipolar diffusion coefficients respectively. The α coefficient controls the cross-over from ambipolar diffusion in the plasma front to free diffusion upstream from the front. It is defined by α = ν i τ M , with ν i the ionization frequency and τ M = ε 0 /[en e (µ e + µ i )] the Maxwell relaxation time.

We checked that propagation speed from this theory is consistent with that observed in the simulations. In order to do this, we needed to average the theoretical propagation speed over the width of the front region because both D ef f and ν i depend strongly on position in the simulations. We defined the beginning of the front as the point x i for which the ambipolar diffusion becomes dominant, that is to say for which ⟨α⟩ y becomes inferior to 1%. The end of the plasma front is therefore located at the position x f for which the ionisation frequency reaches his maximum. The limits of the front are represented by the dashed lines in figure 8. The theoretical front velocity v f is finally defined as :

v f = 1 x f -x i x f xi dx ⟨v(x, y)⟩ y (17)
with v(x, y) = 2 D ef f (x, y)ν i (x, y). The total ionisation frequency ν i is the sum of the three ionization frequencies ν 1 , ν 6 , and ν 7 of reactions 1, 6, and 7, respectively. Equation (17) gives a good approximation for the simulated front velocity as can be seen in figure 9. This velocity depends on the ionization frequency as expected, but also on the diffusion coefficient. Actually the variation of the front velocity is mainly due to the increase of the ionization frequency since the diffusion coefficient does not vary significantly on the range 30 to 40 dBm.

Steady-state results

The plasma-based microwave power limiter was also measured in steady-state configuration, using the experimental setup of figure 3. The power of the microwave signal at the input of the DUT was swept from 25 to 45 dBm, and we measured the input, reflected, and output power ratios. First, it should be noted that without preionization (I DC = 0 mA) no microwave power limitation is observed in simulation or measurement with this input microwave power range. The electric field is not large enough for gas breakdown. This result is consistent with [20]. Hereafter, the MHCD is switched on with a DC current I DC equal to 3 mA.

Figure 10 shows the simulated and measured power ratios in steady-state as a function of the input power. We observe the typical response of a mainly absorptive microwave power limiter. Indeed, all the power ratios remain constant when P in is lower than the limiting threshold of 28 dBm. The limiter is in the "off" state. When P in exceeds 28 dBm, microwave breakdown occurs and the output power ratio decreases while the reflected and dissipated power ratios increase. We note that in the "on" state, the decrease of the output power is mainly due to the absorption of the HPM signal by the plasma. It is consistent with other results obtained for wideband plasma-based microwave power limiter [20].

The simulated and measured power ratios in figure 10 show good overall agreement. Slight discrepancies are observed for low input power. They are mainly due to the fact that the model does not take into account the dielectric losses in the substrate and the coaxial to transmission line transitions which are responsible for a slight impedance mismatch in the DUT. Besides, for input power larger than 40 dBm, simulated and measured power ratios slightly diverge. Two hypotheses can explain these discrepancies. First, during measurements, due to a defect in the assembly of the DUT, the plasma tended to penetrate between the upper part of the dielectric spacer and the dielectric substrate of the suspended microstrip. Such behaviour is not possible in simulation. Secondly, considering that the microwave input power is mainly absorbed by the plasma, thermal effects may occur in the discharge for such high P in values that can affect the local pressure of the gas. Such effects are not taken into account by the model, which may also explain the difference between simulated and measured results at high input power.

The steady state numerical solutions for the ion and electron densities and for the electron temperature are presented in figure 11 as a function of position for the three input powers P in = 30 dBm, P in = 35 dBm and P in = 40 dBm. All these profiles have also been averaged over the plasma domain, along the y axis. The plasma density increases with the incident power as well as the length of the plasma column. The maximum values of ⟨n⟩ y are 3.01 × 10 18 , 1.26 × 10 19 , 3.34 × 10 19 and the length of the plasma columns (the end of the plasma columns are represented by the black dashed lines on figure 11) are of 13.8 mm, 22.8 mm, and 33.6 mm, for each incident power respectively. The electron temperature is uniform in the plasma, and is approximately equal to 3.45 eV, whatever the incident power.

The simulated plasma density is higher than the critical density whatever the incident power. The critical density at 3 GHz is of n c = 1.12 × 10 17 m -3 . The positions for which ⟨n e ⟩ y = n c are represented by the black dashed lines for each profile. These positions correspond approximately to the ones where the quasineutrality is broken, that is to say, to the end of the plasma columns, where the electron temperatures start to decrease.

Conclusion

A plasma-based microwave power limiter has been proposed in a suspended microstrip technology. The incident electric field here ignites a plasma discharge in Argon at 10 Torr that absorbs and reflects the microwave power and thus reduces the level of the signal at the output of the device. The use of a MHCD (Micro Hollow Cathod Discharge) allows to drastically decrease the magnitude of the electric field necessary to ensure microwave gas breakdown. The performances of the limiter have been studied experimentally by measuring the output and absorbed power. These results have also been compared with numerical ones obtained from a self-consistent model that simulates the non linear interaction between the microwave field and the plasma. A simple chemical model have been proposed, that implies four states: neutral argon, effective metastable and resonant states, and one molecular ion. The reaction rates have been calculated with the BOLSIG+ software, with a selfcalculated electron energy distribution function.

The simulations have shown a very good agreement with the experimental results, both in steady state and transient regime, for input powers lower than 40 dBm. These simulations allowed us to accurately link the time evolution of the power ratios with the dynamics of the plasma in the DUT. The propagation of an ionization wave from the MHCD source toward the microwave source is observed both numerically and experimentally. In our plasma discharge, the creation of electrons is mainly due to the step-wise ionization from metastable and resonant states. From numerical results, we deduced that the propagation is due to the gap between the positions of the maximum of the source term for the electrons and the maximum of the electron density. A theoretical model has shown that the propagation velocity depends on the diffusion of the electrons at the head of the ionization front, and on the ionization frequency.

These results allow us to use this model to explore possibles improvement of the system. However, for input powers higher than 40 dBm, the experimental and numerical results diverge and the dynamics of the plasma is not well described by the model. Further studies will be conducted in to understand the causes of this phenomenon.

Figure 1 .

 1 Figure 1. Schematic of the plasma-based microwave power limiter in suspended microstrip technology.

Figure 3 .

 3 Figure 3.Experimental setup for steady-state (blue configuration) and transient (red configuration) high-power microwave measurements.

Figure 4 .

 4 Figure 4. Model geometry of the suspended microstrip circuit. Metal parts are in orange and dielectric parts in grey.

  Simulated electron density.

Figure 5 .

 5 Figure 5. Comparison of simulated and experimental plasma dynamics at different times with P in = 40 dBm at 3 GHz and Argon at 10 Torr.

Figure 6 .

 6 Figure 6. Simulated and measured transient power ratios as a function of time with P in = 40 dBm at 3 GHz and Argon at 10 Torr and I DC = 3 mA.

t 1 t 2 t 3 t 4 n 3 )Figure 7 .

 437 Figure 7. Simulated average electron density as a function of time with P in = 40 dBm at 3 GHz and Argon at 10 Torr.

Figure 8 .

 8 Figure 8. Simulated density profiles and source terms as a function of the position x at t = 75.2 µs, for P in = 40 dBm. The densities and source terms of metastables and resonant states have been summed as well as the two stepwise ionization source terms of electrons.

Figure 9 .Figure 10 .

 910 Figure 9. Simulated, measured and theoretical front velocities as a function of the input power P in at 3 GHz with Argon at 10 Torr and I DC = 3 mA.

PFigure 11 .

 11 Figure 11. Simulated electron and ions density profiles and electron temperature profile as a function of the position for P in ∈ [30.0, 35.0, 40.0] dBm at 3 GHz with Argon at 10 Torr.

Table 1 .

 1 Ar * m,r → e + Ar + Ar + 0.3 × 6.4 × 10 -16 m 3 .s -1 Complete set of reactions taken into account in the numerical model.

	#	Reaction	Reaction rate	Comment	Ref.
	1 2, 3 4, 5 6, 7 8, 9 10 11, 12, 13 Ar * e + Ar → 2e + Ar + e + Ar → e + Ar * m,r e + Ar * m,r → e + Ar e + Ar * m,r → 2e + Ar + e + Ar * m ⇌ e + Ar * r Ar * r → Ar + hν m,r + Ar * m,r → e + Ar + 2 14, 15, 16 Ar * m,r + Ar * pooling ionization EEDF (m 3 .s -1 ) Direct ionization EEDF (m 3 .s -1 ) Impact excitation EEDF (m 3 .s -1 ) De-excitation EEDF (m 3 .s -1 ) Two-step ionization EEDF (m 3 .s -1 ) Metastable quenching Radiative de-excitation 5.85 × 10 5 s -1 Ar * pooling ionization 0.7 × 6.4 × 10 -16 m 3 .s -1 17 Ar + + 2Ar → Ar + 2 + Ar Molecular ion formation 2.5 × 10 -43 m 6 .s -1 18 e + Ar + 2 → Ar + Ar * m 7.35 × 10 -14 × T -2/3 e m 3 .s -1 Dissociative recombination [40] [30] [30] [30] [31-33] [34] [35-37] [38] [38] [39]

Acknowledgments

The authors would like to thank the Délégation Générale pour l'Armement (DGA), the Agence de l'Innovation de Défense (AID), the Région Occitanie Pyrénées-Méditerranée and the Agence Nationale de la Recherche (ANR-21-ASM1-0002) for their financial support.