François Bareille

Julien Wolfersberger

Matteo Zavalloni

Autumn

Institutions and Conservation: the Case of Protected Areas

Keywords: Biodiversity, Democracy, Difference-in-Differences, Natural Parks, Natural Reserves P16, Q56, Q57

We study the link between institutions and the establishment of protected areas for natural resources conservation. Using difference-in-differences we estimate the impacts of democratization on the share of countries' area under protected areas in a panel of 144 countries over 1992-2018. We find that countries that democratized devote higher shares of their land to protected areas, but only in the medium-and long-run. Our preferred estimates indicate that, fifteen years after democratization, the share of the countries' area under PAs increases on average by one percentage point more than what would have been implemented without democratization. We also find evidence that democratization leads countries to implement the different types of protected areas in a non-homogeneous way. Our results are robust to several estimators and treatment definitions, confirming the critical role of institutions for natural resources conservation.

Introduction

While an increasing literature addresses the links between institutions and the environment, both theoretically [START_REF] Fredriksson | Political instability, corruption and policy formation: the case of environmental policy[END_REF][START_REF] Harstad | Conservation contracts and political regimes[END_REF] and empirically [START_REF] Bohn | Ownership risk, investment, and the use of natural resources[END_REF][START_REF] Farzin | Democracy and environmental quality[END_REF][START_REF] Burgess | The political economy of deforestation in the tropics[END_REF], the relationship between institutions and the establishment of protected areas (PAs) has received little attention so far [START_REF] Hu | The political economy of environmental consequences: A review of the empirical literature[END_REF]. This constitutes an important issue since PAs cover about 15% of terrestrial surface and are recognized as a key tool for natural resources conservation [START_REF] Maxwell | Area-based conservation in the twenty-first century[END_REF].

In this paper, we use difference-in-differences to estimate how institutions, proxied by a dichotomous index of democracy, impact the establishment of terrestrial PAs. We unfold two main results related to the size and strictness of established PAs. First, we find that democratizers devote higher shares of land to PAs. Our preferred estimates indicate that, fifteen years after democratization, the share of the countries' area under PAs increases on average by one percentage point (p.p.) more than what would have been established without democratization. Second, we find that democratizers implement the different types of PAs in a non-homogeneous way. Some democratization events are for example associated to a double of the initial share of the countries' area under "wilderness areas" (i.e. large natural reserves), one of the most stringent types of PAs. Others have only significant impacts on the least stringent PA types. Crucially, most of those results are only statistically significant eight to twenty years after democratization, suggesting that institutions affect natural resources conservation mainly in the medium-and long-run.

To reach these results, we build a panel dataset of 144 countries over the 1992-2018 period, which covers the establishment of 68% of all terrestrial PAs. Our primary data source is the World Database on Protected Areas (WDPA) that provides exhaustive information on PAs worldwide.

We measure institutional quality across countries with a democracy index combining Polity5 and Freedom House scores [START_REF] Papaioannou | Democratisation and growth[END_REF][START_REF] Acemoglu | Democracy does cause growth[END_REF]. We then estimate the treatment effects of first democratization event on the size and strictness of PAs using the DID l estimator proposed by De Chaisemartin and D'Haultfoeuille (2022a). In addition to be robust to heterogeneous and dynamic treatment effects, the DID l estimator is particularly adapted to our "unstaggered" setting where countries can switch "in" and "out" democracy. Importantly, we show that our results are robust to alternative estimators, specifications and treatment definitions.

Our analysis contributes to the literature by being the first global-level analysis with an identification strategy that isolates the effect of institutions on PA establishment. Indeed, existing works rely on cross-section regressions [START_REF] Midlarsky | Democracy and the environment: an empirical assessment[END_REF][START_REF] Neumayer | Do democracies exhibit stronger international environmental commitment? a cross-country analysis[END_REF] [START_REF] Van | A semiparametric analysis of determinants of a protected area[END_REF][START_REF] Wurster | Comparing ecological sustainability in autocracies and democracies[END_REF][START_REF] Kashwan | Inequality, democracy, and the environment: A cross-national analysis[END_REF], where several biases prevent causal identification. 1 By comparison, we are using the latest developments in difference-in-differences literature to ensure unbiased measures of the impacts of institutions on PA establishment. Specifically, we are among the first to apply the De Chaisemartin and D'Haultfoeuille (2022a)'s estimator to an applied research question. Given that our estimates are robust to several empirical choices, we believe that they plausibly represent causal impacts of institutions on PA establishment. Our work is also related to a large literature focusing on PA effectiveness (e.g., [START_REF] Joppa | On the protection of protected areas[END_REF], particularly regarding the role of institutions [START_REF] Abman | Rule of law and avoided deforestation from protected areas[END_REF]; Bonilla-Mejía and Higuera-Mendieta, 2019). We view our analysis as complementary to those as we focus on establishment rather than effectiveness. Finally, we contribute to the larger literature on institutions and natural resources conservation where policy implementation has been overlooked [START_REF] Dasgupta | The influence of institutions, governance, and public opinion on the environment: Synthesized findings from applied econometrics studies[END_REF].

The paper is structured as follows. Section 2 presents the data and the econometric strategy.

Section 3 presents the results. Section 4 concludes.

Data and methods

Data description and definitions

Protected Areas. Denoting by A i,t the area under terrestrial PAs in country i and year t and S i the country i's terrestrial size, we define the outcome of interest as the share of a country's area under PAs and write it as a i,t = A i,t /S i . We measure A i,t using the WDPA. This database compiles information on all areas that meet the International Union for the Conservation of Nature (IUCN) definition of a PA, that is "a clearly defined geographical space, recognised, dedicated and managed, through legal or other effective means, to achieve the long-run conservation of nature with associated ecosystem services and cultural values" [START_REF] Dudley | Guidelines for applying protected area management categories[END_REF]. Fed by PA managers, the WDPA is updated every month to provide the latest overview on active PAs [START_REF] Bingham | Sixty years of tracking conservation progress using the world database on protected areas[END_REF]. 2 We use the December 2020 update of the WDPA, which includes 261,004 active PAs. Each observation contains PA localization, establishment date and IUCN classification (within seven categories). This classification categorizes PAs according to their objectives (UNEP-WCMC et al., 1 One exception may be Fouqueray and Papyrakis (2019) using panel data on marine PAs. However, 93.5% of PAs are terrestrial (UNEP-WCMC et al., 2018).

2 The database is available for download at www.protectedplanet.net. Each WDPA update introduces newly designed PAs and removes the recent degazetted sites. These degazetted events could bias our estimates if some of the established PAs after 1992 has been degazetted since. Fortunately, the PADDDtracker.org Version 2.0 database (available at http://padddtracker.org/) indicates that there have been only 85 degazetted PAs and 104 downsized PAs globally between 1992 and 2018, i.e. less than 0.1% of all established PAs. Degazettements are thus unlikely to threat our identification strategy. 2018), from "strict natural reserve" (category Ia) to "PA with sustainable use of natural resources" (category VI). The numbered categories indicate decreasingly stringent constraints on human and economic activities. We call A i,j,t the area under PAs of IUCN category j in country i and year t and a i,j,t = A i,j,t /S i the corresponding share of the country's area under PAs of IUCN category j.

Institutions. Our treatment variable is an indicator measuring democracy. Following the literature (e.g. [START_REF] Acemoglu | Democracy does cause growth[END_REF], we create a dichotomous variable integrating the Revised Combined Polity Score (Polity2 ) from the Polity5 project (Marshall and Gurr, 2020) with the Freedom House's score [START_REF] Repucci | Freedom in the world 2020[END_REF]. The former is a widely used index ranging countries from extreme autocracy (-10) to perfect democracy (+10). It offers the advantage of a yearly categorical variable, available over long time-periods for almost all countries worldwide. The latter categorizes countries in "not free", "partially free" and "free". Specifically, we define a country i in year t as a democracy if Freedom House designates i as "free" or "partially free" and Polity5 assigns a Polity2 score p i,t above p = k. The choice of p is not straightforward. Some authors choose p = 0 (e.g. [START_REF] Acemoglu | Democracy does cause growth[END_REF], while others test greater thresholds (e.g. [START_REF] Papaioannou | Democratisation and growth[END_REF][START_REF] Caselli | Resource windfalls, political regimes, and political stability[END_REF]. Rather than choosing between these alternatives, we systematically define country i as a democracy in year t (D k i,t = 1) if p i,t ≥ k, with k ∈ {0, 2, 4, 6, 8}. We code the country as D k i,t = 0 otherwise. We thus consider five levels of treatment D k in our study (k ∈ {0, 2, 4, 6, 8}), D k i,t being its dichotomous measure for each country-year unit. These measures allow us to test whether the impact of democratization on PA establishment depends on the countries' democratization level.

Panel construction. Both the WDPA and Polity5 database start in 1800 and cover a largethough different -set of countries. Freedom House's score is however only available since 1973. Also, countries and institutions have largely changed since the USSR dissolution. To focus on the contemporary global institutional landscape, we restrict our panel to the 1992-2018 period. About 68% of active terrestrial PAs have been set in this period (see Figure ?? in Appendix). Aggregating information at the country level, we obtain a panel of 144 countries over 27 years. We add measures of countries' land area (km 2 ) and GDP per capita (2017 constant $US) from World Development Indicators database (World Bank, 2020).

Summary statistics. Table 1 provides some summary statistics of our panel. The area under

PAs is 60,511.5 km 2 per country on average over 1992-2018. In 2018, PAs covered 13.1% of countries' terrestrial area (8.3% in 1992). Figure A2 in Appendix shows the location of PAs worldwide at these dates. PAs increased everywhere, including in the Southern hemisphere where the bulk of the remaining primary forests is located. ), and that about 17% of the countries experienced democracy reversals. These "back and forth" events between democracy and non-democracy result in a complex "unstaggered" setting where the "control" and "treated" groups change over time.

We explain in the following how we exploit these differential switches "in" and "out" democracy across countries over time to assess the impacts of institutions on PA establishment.

Econometric strategy

To estimate the impacts of institutions on PA establishment, we rely on a difference-in-differences strategy exploiting differential changes in countries' democratic status (the dichotomous treatment variable D k i,t ) and their share of area under PAs (the outcome variable a i,t ) across countries. We consider two alternative estimation strategies to proceed. Two-way fixed effects. A first strategy is to compare the share of the countries' area under PAs of democratic and non-democratic countries using OLS with two-way fixed effects (TWFE) to control for country-specific and year-specific shocks. The typical TWFE specification to estimate is:

a i,t = α k D k i,t + γ i + λ t + i,t (1)
where a i,t is the share of country i's area under PAs in year t, D k i,t is the dichotomous measure of democratic status, α k is the impact per unit of treatment D k on the share of countries' area under PAs, γ i and λ t are country and year fixed effects respectively and i,t is the remaining error term.

The TWFE estimate αk has been interpreted for long as the average treatment effect (ATE) of a positive value of the dichotomous treatment variable on the outcome variable. While this interpretation is accurate in a simple setting with two groups (control and treated groups) and two periods (pre-and post-treatment periods),3 its validity in more complex designs involving multiple groups and periods is questioned by the recent literature on difference-in-differences (see De Chaisemartin and D'Haultfoeuille, 2022b, for a recent survey). This literature shows that, to interpret αk as an ATE, the setting needs to verify both (i) common trends between control and treated countries and (ii) constant treatment effects across countries and periods (De Chaisemartin and D'Haultfoeuille, 2020; Goodman-Bacon, 2021). 4 While the former has good chances to be valid (and can be formally tested), the latter has been often considered as implausible (Goodman-Bacon, 2021; De Chaisemartin and D'Haultfoeuille, 2022a,b). For this reason, we only estimate equation (1) using OLS in robustness checks (for the five treatment definitions D k with k ∈ {0, 2, 4, 6, 8}) and prefer alternative difference-in-differences estimators for our main analyses.

Difference-in-differences. Given the limits of TWFE to recover treatment effects in complex settings, a second strategy to estimate the impacts of institutions on PA establishment is to turn towards difference-in-differences estimation procedures that are robust to heterogeneous treatment effects. The recent difference-in-differences literature proposes several estimators that are robust to such effects in a multi-period setting (Callaway and Sant'Anna, 2021; Goodman-Bacon, 2021; Sun and Abraham, 2021; De Chaisemartin and D'Haultfoeuille, 2022a). Most of these estimators are however only valid in the context of staggered designs, where the countries received the treatment once and for all (De Chaisemartin and D'Haultfoeuille, 2022b). Only the De Chaisemartin and D'Haultfoeuille (2022a)'s DID l estimator is robust to heterogeneous treatment effects when the design is not staggered. Given our complex setting (see Appendix B), we adopt this approach. Specifically, we aim to measure the dynamic treatment effects of turning in democracy for the first time in year t (D k i,t-1 = 0 and D k i,t = 1) on the share of countries' area under PAs over time a i,t+s . We note β k,s the average treatment effects on the treated (ATT) of a first democratization on PA establishment s years after this event. The idea of De Chaisemartin and D'Haultfoeuille (2022a) to obtain βk,s is to estimate separately the dynamic treatment effects for countries switching "in" and for those switching "out" of treatment for each couple t × s (year of first switch × period after first switch) and to aggregate them afterwards. In the first step, De Chaisemartin and D'Haultfoeuille (2022a) define a difference-in-differences estimator comparing the evolution of the potential outcomes between t -1 and t + s in (i) countries switching of treatment for the first time in t (the "switchers") and (ii) in countries remaining untreated until t (the control groups, including the "never-treated" and "not-yet-switchers"). Under validity of the common trends assumption, 5

De Chaisemartin and D'Haultfoeuille (2022a) note that this estimator is an unbiased measure of the effect of the first treatment switch s periods after it took place in t. In the second step, De Chaisemartin and D'Haultfoeuille (2022a) aggregate these estimators into single coefficients treatment values (e.g. [START_REF] Acemoglu | Democracy does cause growth[END_REF]. In this specification, the coefficient of the contemporaneous treatment is interpreted as the "instantaneous" effect, while the coefficients of the lagged treatments are interpreted as the "dynamic" effects (De Chaisemartin and D'Haultfoeuille, 2022a). However, [START_REF] Sun | Estimating dynamic treatment effects in event studies with heterogeneous treatment effects[END_REF] showed that this specification leads to biased estimates as each period estimate gets contaminated by effects from other periods.

5 De Chaisemartin and D'Haultfoeuille (2022a) propose a procedure to test the parallel trends assumption, relying on the computation of pre-treatment estimates and their joint comparison between switchers and non-switchers before the switchers switch. β k,s by weighting them according to the number of switchers and non-switchers in each couple t × s, multiplying by minus one the estimators associated to switchers "out". De Chaisemartin and D'Haultfoeuille (2022a) show that βk,s is unbiasedly measured since at least some countries are never treated and others are always treated (as in our case, see Appendix B).6

In addition to estimate dynamic treatment effects, De Chaisemartin and D'Haultfoeuille (2022a) propose an additional estimation procedure that provides weights to aggregate all the previous β k,s into a single coefficient β k reflecting the average "total effect" per unit of treatment. Given that the "total effect" refers to the sum of all the instantaneous and dynamic treatment effects, the interpretation of βk is a close approximation to the usual interpretation of αk . More crucially, βk is robust to dynamic and heterogeneous effects. In the following we estimate the β k,s and β k for the five treatment definitions D k (k ∈ {0, 2, 4, 6, 8}).

Results

Size effects

This section presents our results on the impacts of democratization on the share of the country's area under PAs. We first detail the results for the median treatment D 4 -our preferred analysisbefore showing the results with the other treatments and alternative estimators. Preferred analysis. Figure 1 presents the estimates β4,s for s = -13, -12, ..., 19, 20 using the DID l estimator of De Chaisemartin and D'Haultfoeuille (2022a). The solid black line depicts the ATT over time. The black and grey dashed lines plot respectively the 90% and 95% confidence intervals, which are obtained using 100 bootstrap samples with clusters at the country level. The estimates in Figure 1 call for several comments. First, none of the estimates from t -13 to t -1 shows a statistically significant difference between the treated and non-treated countries. This suggests that there is no differential trend in PA establishment before democratization and that, by consequence, the parallel trends assumption holds. This is confirmed by the result of the joint significance test that all the thirteen pre-treatment estimates are null (p = 0.831).

Second, we do not identify any significant impact of the treatment D 4 in the short-run. Indeed, the estimates remain close to zero until four years. While the estimates indicate that democratization starts to positively affect PA establishment after four years, they remain non-significantly This figure plots the dynamic treatment effects of the first observed democratization event on the share of the country's area under PAs (in percentage point). Democracy is defined as a country coded as "free" or "partially free" in Freedom House with a P olity2 score greater than 4 in Polity5. The x-axis represents the time (in years) relative to the year of first switch (t = 0). The vertical dashed line represent the last year before the first switch (t = -1). The solid black line plots the ATT using the De Chaisemartin and D'Haultfoeuille (2022a)'s DID l estimator, which is robust to heterogeneous treatment effects. The black and grey dashed lines represent respectively the 90 and 95 percent confidence interval, obtained by a cluster bootstrapping procedure of 100 random sampling with replacement. The p-value of the joint significance test for the pre-treatment estimates is 0.831, which fails to reject the null hypothesis of no pre-treatment trends.

different from zero at the 10% level. It is only eight years after democratization that democratizers allocate additional 0.3 p.p. of land to PAs compared to what they would have allocated without democratization (significant at the 10% level). Since the share of countries' area under PAs before democratization is on average 11.0 p.p. among the treated countries, the increase of 0.3 p.p.

corresponds to an increase by 3% relatively to the pre-treatment period. The positive impacts of democratization on PA establishment become even greater afterwards. Indeed, Figure 1 indicates that, fifteen years after democratization, democratizers implement on average 1.0 p.p. more PAs than what they would have implemented without democratization (significant at the 5% level).

The impact of democratization on PA establishment increases up to 1.5 p.p. after twenty years (significant at the 10% level).

The long-run nature of the relationship between democratization and PA establishment is not surprising. Indeed, the creation of new PAs itself might take years between the identification phase to the actual establishment. For example, the establishment of the Masoala National Park (Madagascar) took seven years between the identification of the area of interest in 1990, the beginning of the design study in 1993 and its official gazettement in 1997 [START_REF] Kremen | Designing the masoala national park in madagascar based on biological and socioeconomic data[END_REF]. Moreover, such a long-term process is likely to be exacerbated by the turmoil that often characterizes institutional The estimates confirm overall our main results. 7 Indeed, none of the estimates obtained with the different treatments indicates that democracies significantly reduce PAs in the short-or long-run.

Rather, these estimates confirm that democratizers increase by about 1.0 p.p the area under PAs fifteen years after democratization (significantly non-null at the 5% level for all treatments except D 6 ). The long-run effects are only robust at the 10% level (except for D 6 ). Moreover, Table 2 suggests an interesting pattern regarding the effects of democratization as the treatment definition becomes stricter. From D 0 to D 8 , we find that the long-run impacts are more pronounced over time. The effects at twenty years are for example two times larger for D 8 than D 0 . We also find that the significant positive impacts of democratization appear earlier in time as the democratization level threshold increases. While significant positive effects only appear after six, seven and eight years for treatments D 0 , D 2 and D 4 , 8 they become instantaneously positive for D 8 ( β8,0 = 0.059, with an estimated S.D. of 0.020).

TWFE vs. Difference-in-Differences.

As a first -naive -robustness check, we turn to the TWFE estimation of α k in equation ( 1) with OLS for the five treatments D k (k ∈ {0, 2, 4, 6, 8})

and compare them with the estimated total effect β k obtained with the DID l estimator. Table A3 in Appendix presents the obtained estimates. These estimates call for several comments. First, the TWFE estimates are imprecisely measured. Only α8 is significantly positive at the 10% risk level.

By comparison, only β6 is imprecisely estimated with the DID l estimator. Second, the TWFE estimates are very unstable across treatments, ranging from -0.485 to 2.115. By comparison, the DID l estimates are very stable, ranging from 0.362 to 0.780 (and even 0.588 to 0.780 if we only account for the significant estimates). The imprecise and unstable nature of the TWFE estimates is explained by the presence of numerous negative weights for the five treatments (Table A3). While the single significantly non-null TWFE estimate α8 is positive (and would confirm our main result), the presence of negative weights call for its cautious interpretation as negative weights are likely to lead to biased estimates (De Chaisemartin and D'Haultfoeuille, 2020).

Additional robustness checks.

To further test the robustness of our results, we estimate the impact of democratization using alternative difference-in-differences estimators. Specifically, we follow the doubly robust estimation procedure proposed in [START_REF] Acemoglu | Democracy does cause growth[END_REF] to estimate the 7 Despite some pre-treatment estimates are significantly non-null in treatments D 0 and D 2 (see column (2) in Table 2), the joint significance tests proposed in De Chaisemartin and D'Haultfoeuille (2022a) fail to reject the null hypothesis of no pre-treatment trends for all the five treatments D k (k ∈ {0, 2, 4, 6, 8}). This suggests the validity of the parallel trends assumption and the unbiased estimation of the βk,s .

8 These heterogeneous impacts between treatments do not appear in Table 2 as average treatment effects are aggregated over time horizons of five years.

treatment effects of first democratizations (taking place in t) on the share of the countries' area under PAs in t + s (s = -13, ..., 20). 9 The doubly robust estimation procedure only compares information on democratizers (switchers "in") with information on not-yet and never treated countries. The doubly robust estimator thus do not fully exploit the complex "unstaggered" setting of our panel and, by consequence, the obtained estimates cannot be quantitatively compared to the DID l estimates. The doubly robust estimator has however the advantage to allow the inclusion of time-invariant covariates [START_REF] Acemoglu | Democracy does cause growth[END_REF], which the DID l estimator cannot (De Chaisemartin and D'Haultfoeuille, 2022a). Such additional information is particularly useful to control for unobserved heterogeneity. Specifically, a possible omitted variable that could affect our estimates is the level of economic development as it may be correlated to both institutions [START_REF] Acemoglu | Democracy does cause growth[END_REF] and PAs [START_REF] Van | A semiparametric analysis of determinants of a protected area[END_REF]. For the five levels of treatments D k (k ∈ {0, 2, 4, 6, 8}), we use the doubly robust estimator to estimate the dynamic treatment effects of democratization on the share of treated countries' area under PAs for two models: one including only the lagged outcome variable a i,t-1 and year fixed effects (in the spirit of [START_REF] Acemoglu | Democracy does cause growth[END_REF], the other additionally including the lagged GDP per capita (i.e. before the country democratized). 10 Tables A4 and A5 in Appendix present the estimated ATT for these two models. We find qualitatively similar results to our preferred analysis. Indeed, none of the doubly robust estimates obtained with the different treatments indicates that democratizers significantly reduce PAs in the short-or long-run. The two set of estimates show a similar picture to those in Table 2 using the DID l estimator: greater democratization events are associated with both earlier appearance of significant positive effects and larger long-run effects. 11 In particular, the positive impacts of 9 As explained in [START_REF] Acemoglu | Democracy does cause growth[END_REF], such a semi-parametric approach uses additional covariates to build a doubly robust estimator that combines (i) linear estimates of the potential outcomes of the democratizers (in the spirit of a regression adjustment approach à la Kline, 2011) reweighted by (ii) the scores obtained from a probit model explaining the likelihood to democratize (in the spirit of a propensity score reweighting approach à la [START_REF] Angrist | Semiparametric estimates of monetary policy effects: string theory revisited[END_REF]. The procedure then measures the effect of a first democratization in t on PA establishment in t + s as the difference between (i) the evolution of the reweighted estimated potential outcomes of democratizers between t and t + s and (ii) the evolution of the potential outcomes for the not-yet and never treated countries on a similar period. The doubly robust estimator has the advantage to be consistent if either the regression adjustment or propensity score reweighting approaches is valid [START_REF] Imbens | Recent developments in the econometrics of program evaluation[END_REF]. As the doubly robust estimation procedure is independently implemented for each s (s = -13, .., 20), the obtained estimates are robust to heterogeneous effects.

10 Indeed, the contemporaneous GDP per capita cannot be used directly to condition our estimations as it presents "bad control" properties that could bias the estimation of the dynamic treatment effects.

11 The magnitude of the impacts are larger in Tables A4 and A5 than those identified in Table 2. This is not surprising as the doubly robust estimation procedure presented in [START_REF] Acemoglu | Democracy does cause growth[END_REF] does not exploit the diversity of democratizations. It only exploits information related to first democratization events, ignoring additional information from democratization reversals (and information from countries that remained always in democracy).

democratization on PA establishment are conserved even after controlling for the countries' level of development (Table A5).12

Strictness effects

In addition to the estimation of the impacts of democratization on the share of the country's area under PAs, we also seek to identify potential democratization effects on the types of established PAs (i.e. the different IUCN categories). We thus estimate the impacts of a first democratization D k (k ∈ {0, 2, 4, 6, 8}) occurring in t on the share of the countries' area under PAs in category IUCN j a i,j,t (with j ∈ {Ia, Ib, II, III, IV, V, VI}) using the DID l estimator. Table 3 reports the estimated total effect when the joint significance test on pre-treatment estimates supports the common trends assumption, as well as the average shares of PAs under different IUCN categories the year before the countries democratized (a i,j,t-1 ).

Table 3 shows that no democratization is associated with any decrease in PAs of a particular IUCN category. When precisely estimated, the total effects rather indicate positive impacts of democratization on the shares of PAs under the different IUCN categories. As for the size effects, all these effects grow over time with the largest effects in the long run (see Table A6 in Appendix for an illustration of the dynamic treatment effects for treatment D 0 ). The effects of similar democratizations are however different for the different IUCN categories (Table 3 ). For example, the average total effects of an additional unit of treatment D 0 ranges from 0.011 for the IUCN III category ("Natural monument or feature") to 0.073 for the IUCN V category ("Protected Landscape"). However, the initial shares of PAs under the different categories were different before democratization. When reported to the initial shares, the highest total effect of democratization D 0 actually relates to IUCN category Ib ("Wilderness Area"), one of the strictest category of PAs. 13 The smallest impacts relate to the IUCN category II ("Natural Parks"). Interestingly, the positive impacts on the strictest PAs disappear for stricter democratization definitions.

On the contrary, the largest impacts -both in absolute and relative terms -for D 8 relates to IUCN category VI ("PA with sustainable use of natural resources"), which represents the least stringent PAs. Note. The first part of the table indicates estimates of the total effect of a democratization on the share of the country's area under PA of different IUCN categories, indicated in the column labels. Democracy is defined as a country coded as "free" or "partially free" in Freedom House with a P olity2 score greater than k (k ∈ {0, 2, 4, 6, 8}) in Polity5. We report average total estimates of the ATT using the De Chaisemartin and D'Haultfoeuille (2022a)'s DID l average estimator. Below each estimate we report standard errors (clustered at the country level) obtained via bootstrapping. * and ** indicate a p-value lower than 0.1 and 0.05 respectively. The reported estimates are those for which the joint significance test fails to reject the null hypothesis of no-pre-treatment trends. The symbol "-" indicates the rejection of the parallel trends assumption. The second part of the table indicates the initial share of PAs in the different IUCN categories before countries switch for the first time the treatment D k (k ∈ {0, 2, 4, 6, 8}).

Conclusion

In this paper, we investigate how institutions affect natural resources conservation by looking at terrestrial PA establishment. For this purpose, we built a panel dataset over 27 years and 144 countries (covering 68% of active terrestrial PAs) and exploited the occurrence of democratization events across countries to run difference-in-differences analyses with estimators that are robust to heterogeneous and dynamic treatment effects. We identified two main results. First, we found that democratizers allocate significantly larger share of their land area to PAs. Second, we found that democratizers implement the different types of PAs in a non-homogeneous way. Importantly, these effects mainly appear in the medium-and long-run, i.e. eight to twenty years after democratization.

These results are consistent with historical examples worldwide. Eventually, we also found that higher democratization scores are associated to (i) earlier positive impacts of institutions on PA establishment over time, (ii) larger long-run impacts but (iii) slightly less stringent established PAs.

Because these results are robust to several estimators, specifications and treatment definitions, we believe that they represent plausibly causal impacts of institutions, as measured by democratization, on PA establishment. The question of the mechanisms that conduct democracies to implement larger and more stringent PAs was left out of our analysis. There could be several reasons that explain such changes, such as elections' outcomes, the presence of NGOs in a country or the signature of bilateral or international agreement [START_REF] Groom | Carbon emissions reductions from indonesia's moratorium on forest concessions are cost-effective yet contribute little to paris pledges[END_REF]. We believe that the identification of the role of these underlying mechanisms constitutes an interesting avenue to pursue in future research.

B Democratization events over time

Democratization events are not similar to each others. First, democratization events can differ in their amplitudes. For example, Ethiopia switched in democracy in 1993 according to treatment definition D 0 but remained in non-democracy according to other treatment definitions. Second, democratization events can occur at different periods. For example, Ghana democratized in 1996 and Tanzania democratized in 2015 (using definition D 0 ). Third, democracy can be reversed. On the one hand for example, Burundi democratized in 2003 after a period of eleven years in nondemocracy (using definition D 0 ), but experienced a reversal towards non-democracy in 2014. On the other hand, Burkina Faso democratized in 2001 and has remained a democracy ever since (using definition D 0 ). Tables A1 andA2 summarize this diversity of the democratization events in our panel for the five treatment definitions D k (k ∈ {0, 2, 4, 6, 8}), depending on whether countries switch "in" or "out" democracy.

Regarding Table A1, countries turn more often to democracies than to non-democracies for all treatment definitions. Both democratization and autocratization events are more frequent in the beginning of the panel. They are also about two times more frequent for D 0 , D 2 , D 4 and D 6 than for D 8 in each decade. Overall, out of the 144 countries that constitute our panel, 26 to 41 countries turned to democracies and 16 to 30 reversed, depending on the treatment definition (Table A2). Looking for example to countries switching "in" treatment D 0 (D 0 i,t-1 = 0 and D 0 i,t = 1), 39 countries democratized between 1992 and 2018, 64 remained always democracies and 41 never experienced democratization. This last category includes 31 countries that remained always non-democracies and 10 others that definitely reversed to non-democracies between 1992 and 2018. From the opposite perspective, looking at countries switching "out" D 0 (D 0 i,t-1 = 1 and D 0 i,t = 0), 30 countries experienced democracy reversals (including the 10 countries mentioned above), while 31 countries remained always non-democracies, and 83 countries never experienced democracy reversals (including the 64 mentioned countries remaining democratic between 1992 and 2018). Note. This table indicates estimates of the effect of a democratization on the share of the country's area under PAs over 5-years time horizons, indicated in the column labels. Democracy is defined as a country coded as "free" or "partially free" in Freedom House with a P olity2 score greater than k (k ∈ {0, 2, 4, 6, 8}) in Polity5. We report doubly robust estimates of the ATT using the procedure of Acemoglu et al. (2019), using lagged share of countries' area under PAs, lagged GDP per capita and year fixed effects as covariates.

Below each estimate we report standard errors (clustered at the country level) obtained via bootstrapping. *, **, *** indicate a p-value lower than 0.1, 0.05 and 0.01 respectively. We report the p-value of the joint significance tests on the pre-treatment estimates between hooks. Note. This table indicates estimates of the effect of a first democratization on the share of the country's area under PAs over 5-years time horizons, indicated in the column labels. Democracy is defined as a country coded as "free" or "partially free" in Freedom House with a P olity2 score greater than 0 in Polity5. We report estimates of the ATT using the De Chaisemartin and D'Haultfoeuille (2022a)'s DID l estimator. Below each estimate we report standard errors (clustered at the country level) obtained via bootstrapping. * and ** indicate a pvalue lower than 0.1 and 0.05 respectively. We report the p-value of the joint significance tests on the pre-treatment estimates between hooks.

E Strictness effects: dynamic treatment effects of a first democratization event on the share of the countries' area under PAs of different IUCN categories

Figure 1 :

 1 Figure 1: Estimates of the effects of democratization on the share of the country's area under PAs.

  changes. For example, Madagascar officially democratized in 1991 after socialist regime's dissolution. In 2003, twelve years after democratization, the Malagasy government committed to triple the area under PAs in the country by 2008 (Gardner et al., 2018). Similar long-term patterns have been identified in the 1980s in Portugal and Guatemala, where the establishment of first PAs took place five to ten years after the first elections (Sundberg, 2003; Pinto and Partidário, 2012).Democratization score. The impacts of democratization on PA establishment may depend on its definition. Table2displays the dynamic treatment effects obtained with the DID l estimator for the five treatment definitions and aggregated on successive periods of five years.

Figure A1 :

 A1 Figure A1: Cumulative number and area under protected areas in time worldwide (Source: WDPA and authors' computations). The black lines indicate the evolution of the number and cumulative area under PAs from 1900 to 2018. The vertical red lines indicate the beginning of our panel, in 1992.

Figure A2 :

 A2 Figure A2: Evolution of the share of countries' area under PAs worldwide between 1992 and 2018 (Source: WPDA Project)

Figure A3 :

 A3 Figure A3: Evolution of Polity2 scores worldwide between 1992 and 2018 (Source: Polity5 Project)

Table 1 :

 1 Table 1 further indicates that the most common PAs are National Parks (IUCN category II), which occupy about 31.6% of all PAs. The least frequent PAs are "natural monuments" (IUCN category III), which occupy 1.5% of all PAs. Summary statistics (N=3,662)

	Variables	Mean	S.D.	Min	Max
	Cumulative terrestrial area under PAs (km 2 ) 60,511.5	154,028.4	0.3	1,250,856.4
	IUCN category Ia (km 2 )	4,532.7	23,543.6	0.0	228,033.2
	IUCN category Ib (km 2 )	6,270.0	39,229.8	0.0	390,237.9
	IUCN category II (km 2 )	15,763.8	40,806.5	0.0	353,470.2
	IUCN category III (km 2 )	1,088.3	6,148.5	0.0	63.656.46
	IUCN category IV (km 2 )	11,390.6	51,846.9	0.0	426,102.2
	IUCN category V (km 2 )	6,008.3	25,308.0	0.0	251,233.0
	IUCN category VI (km 2 )	12,451.8	47,755.7	0.0	365,396.5
	Non-attributed IUCN category (km 2 )	3,005.8	9,604.1	0.0	105,982.4
	Additional area under PAs per year (km 2 )	438.0	1,643.0	0.0	32,014.2
	Country terrestrial size (km 2 )	832,572.2 2,083,846.0 670.3 16,389,950.1
	Revised Combined Polity Score (P olity2)	3.6	6.4	-10.0	10.0
	GDP per capita (1,000 constant 2017$)	12.1	17.9	0.2	112.0
	Regarding institutions, Table 1 indicates that the average P olity2 score is 3.6 over 1992-2018,
	underlying that most countries are democracies. Precisely, the average P olity2 score increased
	from 2.2 in 1992 to 4.2 in 2018, with a peak at 4.4 in 2015. Appendix B provides elements on the
	diversity of democratization events in our panel for each treatment definition D k (k ∈ {0, 2, 4, 6, 8}).
	Specifically, Figure A3 shows the P olity2 score in 1992 and 2018 across countries, suggesting that
	democratizations mostly occurred in Africa. Table A1 indicates the number of democratizations
	and reversals per decade. It shows that, in our panel, there have been on average about two de-
	mocratizations per year for treatments D 0 to D 6 but only one D 8 democratization per year. It also
	indicates that democratizations have been more frequent in the nineties. Table A2 displays the
	number of countries switching respectively to democracies and non-democracies according to treat-
	ment definitions. It indicates that about 25% of the countries experienced democratization during
	the 1992-2018 period (except for treatment D 8

Table 2 :

 2 Estimates of the dynamic treatment effects of first democratization event on the share of the countries' area under PAs depending on the threshold on P olity2 score Note. This table indicates estimates of the effect of a democratization on the share of the country's area under PAs over 5-years time horizons, indicated in the column labels. Democracy is defined as a country coded as "free" or "partially free" in Freedom House with a P olity2 score greater than k (k ∈ {0, 2, 4, 6, 8}) in Polity5. We report estimates of the ATT using the De Chaisemartin and D'Haultfoeuille (2022a)'s DID l estimator. Below each estimate we report standard errors (clustered at the country level) obtained via bootstrapping.

				Average Effects from			Joint test
		-13 to -10	-9 to -5	-4 to 0	1 to 5	6 to 10	11 to 15 16 to 20	for
		Years	Years	Years	Years	Years	Years	Years	Parallel trends
		(1)	(2)	(3)	(4)	(5)	(6)	(7)	Assumption
			A. Treatment D 0 ("free" or "partially free" and P olity2 ≥ 0)
	Avg. Effect on	-0.176	-0.427 ** -0.097	0.354	0.734 ** 1.018 ** 0.978 **	Valid
	share of PAs	(0.447)	(0.206)	(0.069) (0.226)	(0.391)	(0.479)	(0.480)	[p=0.554]
			B. Treatment D 2 ("free" or "partially free" and P olity2 ≥ 2)
	Avg. Effect on	-0.196	-0.471 ** -0.111	0.364	0.731 *	1.000 ** 0.806 *	Valid
	share of PAs	(0.469)	(0.225)	(0.074) (0.222)	(0.412)	(0.510)	(0.488)	[p=0.350]
			C. Treatment D 4 ("free" or "partially free" and P olity2 ≥ 4)
	Avg. Effect on	-0.098	-0.317	-0.083	0.036	0.255 *	0.728 ** 1.209 *	Valid
	share of PAs	(0.407)	(0.268)	(0.062) (0.070)	(0.149)	(0.367)	(0.655)	[p=0.831]
			D. Treatment D 6 ("free" or "partially free" and P olity2 ≥ 6)
	Avg. Effect on	-0.465	-0.285	-0.042	0.025	0.102	0.270	0.809	Valid
	share of PAs	(0.404)	(0.172)	(0.055) (0.082)	(0.194)	(0.484)	(1.055)	[p=0.746]
			E. Treatment D 8 ("free" or "partially free" and P olity2 ≥ 8)
	Avg. Effect on	0.016	-0.205	-0.061 0.196 ** 0.598 *** 1.071 ** 1.794 *	Valid
	share of PAs	(0.786)	(0.269)	(0.054) (0.082)	(0.280)	(0.510)	(0.971)	[p=0.789]

*, **, *** indicate a p-value lower than 0.1, 0.05 and 0.01 respectively. We report the p-value of the join significance tests on the pre-treatment estimates between hooks.

Table 3 :

 3 Estimates of the average total effect of democratization on the share of the countries' area under PAs of different IUCN categories depending on the threshold on P olity2 score

				Average Total Effect For
		IUCN Ia IUCN Ib IUCN II IUCN III IUCN IV IUCN V	IUCN VI
		(1)	(2)	(3)	(4)	(5)	(6)	(7)
		A. Treatment D 0 ("free" or "partially free" and P olity2 ≥ 0)
	Avg. Effect on	0.013 ** 0.024 ** 0.072 *	0.011 **	0.037 *	0.073 **	0.047 *
	share of PAs	(0.006)	(0.010)	(0.042)	(0.005)	(0.022)	(0.037)	(0.024)
		B. Treatment D 2 ("free" or "partially free" and P olity2 ≥ 2)
	Avg. Effect on	0.020 ** 0.021 **	0.064	0.010 **	-	0.062	0.053 **
	share of PAs	(0.009)	(0.009)	(0.040)	(0.004)		(0.044)	(0.025)
		C. Treatment D 4 ("free" or "partially free" and P olity2 ≥ 4)
	Avg. Effect on	0.013	0.031	0.069	0.005	-0.057	0.097	0.027
	share of PAs	(0.010)	(0.022)	(0.052)	(0.010)	(0.041)	(0.059)	(0.033)
		D. Treatment D 6 ("free" or "partially free" and P olity2 ≥ 6)
	Avg. Effect on	0.018	0.019	0.022	-	-	0.019	-
	share of PAs	(0.018)	(0.029)	(0.056)			(0.048)
		E. Treatment D 8 ("free" or "partially free" and P olity2 ≥ 8)
	Avg. Effect on	-	0.017	-	0.020 **	0.060	0.070	0.052 *
	share of PAs		(0.010)		(0.010)	(0.038)	(0.044)	(0.030)
				Average Initial Share For
		IUCN Ia IUCN Ib IUCN II IUCN III IUCN IV IUCN V	IUCN VI
	First switch "in" D 0	0.342	0.019	3.059	0.249	1.467	0.179	1.701
	First switch "in" D 2	0.390	0.017	2.827	0.248	1.342	0.188	1.460
	First switch "in" D 4	0.533	0.310	2.886	0.278	1.375	0.179	1.386
	First switch "in" D 6	0.531	0.617	2.952	0.124	1.732	0.103	0.613
	First switch "in" D 8	0.417	0.674	1.991	0.282	0.905	0.136	0.649

Table A1 :

 A1 Annual average of the number of democratization and autocratization events per decades (Source: authors' own computations based on Polity5 and Freedom House databases)

			Year of events (annual average per decades)
			Switching "in"			Switching "out"
		D 0	D 2	D 4	D 6	D 8	D 0	D 2	D 4	D 6	D 8
	90's	2.43 2.29 2.14 1.86 1.14	1.71 1.57 1.86 1.00 1.00
	00's	1.82 1.91 2.09 2.00 1.36	1.40 1.30 1.30 1.30 0.60
	10's	1.44 1.67 1.67 1.56 0.56	1.22 1.22 1.22 1.11 0.44
	Average 1.81 1.92 2.00 1.85 1.04	1.42 1.35 1.42 1.15 0.65
	Note. The Table indicates the average number of democratization events (D k i,t-1
	= 0 and D k i,t = 1) and autocratization events (D k i,t-1 = 1 and D k i,t = 0) per year for
	each decade of our panel and each treatment definition.

Table A2 :

 A2 Number of switchers "in" and "out" per treatment during the 1992-2018 period (Source: authors' own computations based on Polity5 and Freedom House databases)

					Group Size per Treatment
			Switching "in"			Switching "out"
		D 0 D 2 D 4 D 6 D 8	D 0 D 2 D 4 D 6	D 8
	Never-treated	41	44	46	54	79	83	69	79	77	60
	Treated	39	40	41	39	26	30	26	28	24	16
	Always-treated 64	60	57	51	39	31	35	37	43	68
	Total	144 144 144 144 144	144 144 144 144	144
	Note. The table displays the number of countries that get never-, always-or some-
	times treated according to treatment definitions D k (k ∈ {0, 2, 4, 6, 8}). Table displays
	in separated columns the countries switching "in" democracy (D k i,t-1 = 0 and D k i,t = 1)
	and those switching "out" democracy (D k i,t-1 = 1 and D k i,t = 0).

Table A3 :

 A3 Average estimate of the effect of democratization on the share of PAs using TWFE (OLS) and DID l This table indicates average estimates of the effect of democratization on the share of the country's area under PA. Democracy is defined as a country coded as "free" or "partially free" in Freedom House with a P olity2 score greater than k (k ∈ {0, 2, 4, 6, 8}) in Polity5. We report average estimates obtained with OLS and DID l . Below each estimate we report standard errors (clustered at the country level), obtained via bootstrapping for the DID l estimates. *, **, *** indicate a p-value lower than 0.1, 0.05 and 0.01 respectively. Below each standard errors we report the number of ATTs receiving negative weights in the TWFE estimation procedure, expressed out of the total number of ATTs used in the TWFE estimation (obtained with the twowayfeweights Stata command from De Chaisemartin and D'Haultfoeuille, 2020).

			Average effect
		Two-Way Fixed Effects	Difference-in-differences (DID l )
		A. Treatment D 0 ("free" or "partially free" and P olity2 ≥ 0)
	Avg. Effect on	-0.491	0.780 **
	share of PAs	(0.449)	(0.378)
	ATTs with neg. weights	983/2462
		B. Treatment D 2 ("free" or "partially free" and P olity2 ≥ 2)
	Avg. Effect on	-0.427	0.742 *
	share of PAs	(0.465)	(0.379)
	ATTs with neg. weights	869/2392
		C. Treatment D 4 ("free" or "partially free" andP olity2 ≥ 4)
	Avg. Effect on	-0.485	0.588 *
	share of PAs	(0.451)	(0.306)
	ATTs with neg. weights	934/2316
		D. Treatment D 6 ("free" or "partially free" and P olity2 ≥ 6)
	Avg. Effect on	0.386	0.362
	share of PAs	(0.591)	(0.600)
	ATTs with neg. weights	794/2035
		E. Treatment D 8 ("free" or "partially free" and P olity2 ≥ 8)
	Avg. Effect on	2.115 *	0.765 **
	share of PAs	(1.250)	(0.335)
	ATTs with neg. weights	671/1496
	Note.

D Robustness

check: doubly-robust estimator using the Ace- moglu et al. (2019)'s procedureTable A4 :

 A4 Doubly robust estimates of the effects of democratization on the share of the country's area under PAs depending on the threshold on P olity2 score (conditions on lagged share of countries' area under PAs and year fixed effects)Note. This table indicates estimates of the effect of a democratization on the share of the country's area under PAs over 5-years time horizons, indicated in the column labels. Democracy is defined as a country coded as "free" or "partially free" in Freedom House with a P olity2 score greater than k (k ∈ {0, 2, 4, 6, 8}) in Polity5. We report doubly robust estimates of the ATT using the procedure of[START_REF] Acemoglu | Democracy does cause growth[END_REF], using lagged share of countries' area under PAs and year fixed effects as covariates. Below each estimate we report standard errors (clustered at the country level) obtained via bootstrapping. *, **, *** indicate a p-value lower than 0.1, 0.05 and 0.01 respectively. We report the p-value of the joint significance tests on the pre-treatment estimates between hooks.

				Average Effects from			Joint test
		-13 to -10 -9 to -5 -4 to 0	1 to 5	6 to 10 11 to 15 16 to 20 for Conditional
		Years	Years	Years	Years	Years	Years	Years	Parallel trends
		(1)	(2)	(3)	(4)	(5)	(6)	(7)	Assumption
		A. Treatment D 0 ("free" or "partially free" and P olity2 ≥ 0)
	Avg. Effect on	-0.073	0.052	0.006	0.048	0.530	2.045	3.170	Valid
	hare of PAs	(0.285)	(0.108) (0.043) (0.102) (0.482)	(1.354)	(2.381)	[p=0.398]
		B. Treatment D 2 ("free" or "partially free" and P olity2 ≥ 2)
	Avg. Effect on	-0.139	-0.01	-0.049	0.088	0.693	2.772 *	4.042 *	Valid
	hare of PAs	(0.264)	(0.086) (0.041) (0.118) (0.578)	(1.631)	(2.354)	[p=0.421]
		C. Treatment D 4 ("free" or "partially free" and P olity2 ≥ 4)
	Avg. Effect on	-0.357	-0.256	-0.084	0.101	0.707	2.593 ** 3.927 *	Unvalid
	share of PAs	(0.408)	(0.255) (0.073) (0.131) (0.562)	(1.263)	(2.069)	[p=0.065]
		D. Treatment D 6 ("free" or "partially free" and P olity2 ≥ 6)
	Avg. Effect on	-0.819	-0.327	-0.086	0.248	1.317	2.219	0.406	Unvalid
	share of PAs	(0.542)	(0.269) (0.065) (0.232) (0.989)	(1.398)	(1.995)	[p=0.098]
		E. Treatment D 8 ("free" or "partially free" and P olity2 ≥ 8)
	Avg. Effect on 0.304 ***	-0.130	0.033	1.336 * 3.632 *	4.008	7.846	Valid
	share of PAs	(0.116)	(0.365) (0.142) (0.800) (2.173)	(2.788)	(5.971)	[p=0.102]

Table A5 :

 A5 Doubly robust estimates of the effects of democratization on the share of the country's area under PAs depending on the threshold on P olity2 score (conditions on lagged share of countries' area under PAs, lagged GDP per capita and year fixed effects)

				Average Effects from			Joint test
		-13 to -10 -9 to -5 -4 to 0	1 to 5 6 to 10 11 to 15 16 to 20	for Conditional
		Years	Years	Years	Years	Years	Years	Years	Parallel trends
		(1)	(2)	(3)	(4)	(5)	(6)	(7)	Assumption
			A. Treatment D 0 ("free" or "partially free" and P olity2 ≥ 0)
	Avg. Effect on	-0.113	0.020	0.001	0.069	0.594	2.328	4.177	Valid
	share of PAs	(0.261)	(0.08)	(0.044) (0.102) (0.499) (1.501)	(2.788)	[p=0.404]
			B. Treatment D 2 ("free" or "partially free" and P olity2 ≥ 2)
	Avg. Effect on	-0.174	-0.047	-0.056	0.102	0.746	2.838 *	4.066 *	Valid
	share of PAs	(0.239)	(0.077) (0.042) (0.118) (0.571) (1.575)	(2.238)	[p=0.405]
			C. Treatment D 4 ("free" or "partially free" and P olity2 ≥ 4)
	Avg. Effect on	-0.391	-0.301	-0.09	0.104	0.705	2.545 *	3.885 *	Valid
	share of PAs	(0.406)	(0.257) (0.064) (0.138) (0.627) (1.484)	(2.280)	[p=0.346]
			D. Treatment D 6 ("free" or "partially free" and P olity2 ≥ 6)
	Avg. Effect on	-0.824 *	-0.371 * -0.094	0.255	1.321	2.133	0.073	Unvalid
	share of PAs	(0.464)	(0.214) (0.063) (0.219) (0.875) (1.368)	(2.233)	[p=0.074]
	Avg. Effect on	0.270	-0.138	0.039	1.299	3.511 *	3.825	7.355	Valid
	share of PAs	(0.162)	(0.395) (0.162) (0.883) (2.036) (2.610)	(4.936)	[p=0.122]

E. Treatment D 8 ("free" or "partially free" and P olity2 ≥ 8)

Table A6 :

 A6 Estimates of the dynamic treatment effects of a first democratization event on the share of the countries' area under PAs of different IUCN categories

				Average Effects from			Joint test
		-13 to -10 -9 to -5 -4 to 0	1 to 5	6 to 10 11 to 15 16 to 20	for
		Years	Years	Years	Years	Years	Years	Years	Parallel trends
		(1)	(2)	(3)	(4)	(5)	(6)	(7)	Assumption
				A. IUCN Ia ("Strict Nature Reserve")
	Avg. Effect on	0.254	0.002	-0.003	0.006	0.013 ** 0.019 ** 0.018 *	Valid
	share of PAs	(0.201)	(0.045) (0.002) (0.004)	(0.006)	(0.008)	(0.011)	[p=0.431]
				B. IUCN Ib ("Wilderness Area")
	Avg. Effect on -0.031 *	-0.024	-0.006 0.011 ** 0.025 ** 0.032 ** 0.034 **	Valid
	share of PAs	(0.016)	(0.015) (0.004) (0.005)	(0.012)	(0.015)	(0.015)	[p=0.437]
				C. IUCN II ("Natural Park")
	Avg. Effect on	0.103	-0.066	-0.020	0.037	0.078	0.091 *	0.088	Valid
	share of PAs	(0.125)	(0.052) (0.017) (0.028)	(0.048)	(0.051)	(0.058)	[p=0.377]
			D. IUCN III ("Natural Monument or Feature")
	Avg. Effect on	-0.007	-0.006	-0.002	0.005	0.011 ** 0.016 ** 0.017 *	Valid
	share of PAs	(0.011)	(0.006) (0.002) (0.003)	(0.005)	(0.007)	(0.009)	[p=0.109]
			E. IUCN IV ("Habitat or Species Management Area")
	Avg. Effect on	-0.010	-0.019	-0.010	0.018	0.042 *	0.043	0.043	Valid
	share of PAs	(0.039)	(0.022) (0.008) (0.014)	(0.024)	(0.036)	(0.043)	[p=0.129]
			F. IUCN V ("Protected Landscape or Seascape")
	Avg. Effect on	-0.087	-0.088 * -0.019	0.039	0.076 *	0.092	0.103	Valid
	share of PAs	(0.058)	(0.049) (0.015) (0.025)	(0.046)	(0.061)	(0.068)	[p=0.232]
				G. IUCN VI ("PA with Sustainable Use of Natural Resources")
	Avg. Effect on	-0.022	-0.020	-0.004	0.019	0.046 *	0.069 *	0.070 *	Valid
	share of PAs	(0.081)	(0.027) (0.005) (0.013)	(0.027)	(0.041)	(0.04)	[p=0.789]

This interpretation of TWFE estimates is valid since the two groups would have evolved in parallel without the treatment, i.e. under the verification of the "parallel trends" assumption.

There are two main possible source of heterogeneity in the treatment effects: across countries and across periods. Regarding the former, equation (1) implicitly assumes that the treatment effects are identical across countries within a similar group (control or treated). However, this assumption is problematic as the control group in equation (1) includes both the never-and always-treated countries. De Chaisemartin and D'Haultfoeuille (2020) show that, given TWFE estimates are weighted sums of treatment effects in each country-year unit, some outcome variable changes can be attributed negative weights for always-treated countries. In case of heterogeneous treatment effects, negative weights lead to biased TWFE estimates (De Chaisemartin and D'Haultfoeuille, 2020).Regarding the latter, equation (1) assumes that countries' potential outcomes only depend on their contemporaneous treatment, not on their past treatments (Goodman-Bacon, 2021). The shape of the impacts of the treatment can be very heterogeneous over time yet: cumulative, concentrated in the years following the first treatment change or, instead, appearing after a sufficiently long period. They can even be positive in the short run but negative in the long run (or vice versa). The usual solution to estimate such dynamic effects in a TWFE setting is to add lagged

We determine the number of pre-and post-treatment periods accordingly (13 and 20 respectively).

The additional condition on the development level actually removes some pre-trends between treated and nontreated countries, improving the validity of the conditional parallel trends assumption.

 13 IUCN defines PAs of category Ib as an "area generally larger and protected in a slightly less stringent manner than that of strict nature reserves".


Emilien Prost
email: emilien.prost@univ-lyon1.fr

Vincent Loriot

Eric Constant

Isabelle Compagnon

Luc Bergé

Franck Lépine

Stefan Skupin

Air-photonics based terahertz source and detection system

Keywords: terahertz generation, plasma, filamentation, ultrashort pulse propagation, spectroscopy

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Terahertz (THz) radiation lies in the part of the electromagnetic spectrum usually defined between 0.1 THz and 30 THz. As THz photons are weakly energetic, they do not ionize matter and can only excite free charges, optical phonons in solids and molecular vibrations and rotations occurring over ∼ ps time scales. The past decades have known tremendous technological developments in the elaboration of THz emitters [START_REF] Tonouchi | Cutting-edge terahertz technology[END_REF][START_REF] Redo-Sanchez | Review of terahertz technology readiness assessment and applications[END_REF][START_REF] Zhang | Intense terahertz radiation: generation and application[END_REF], particularly linked to the advent of ultrafast laser sources. This technological progress opened up the route of new diagnostics for medical imaging and cancer detection [START_REF] Yu | The potential of terahertz imaging for cancer diagnosis: A review of investigations to date[END_REF], new devices for molecular spectroscopy [START_REF] Bergé | Terahertz spectroscopy from air plasmas created by two-color femtosecond laser pulses: The ALTESSE project[END_REF] as well as for greenhouse gas detection and plant life monitoring [START_REF] Gente | Monitoring leaf water content with THz and sub-THz waves[END_REF]. More recently, the production of THz waves operating at high field strengths attracted interest for various applications, e.g., in solid-state physics (see, e.g., the recent review [START_REF] Salén | Matter manipulation with extreme terahertz light: Progress in the enabling thz technology[END_REF]), to probe water dynamics [START_REF] Novelli | Towards Intense THz Spectroscopy on Water: Characterization of Optical Rectification by GaP, OH1, and DSTMS at OPA Wavelengths[END_REF], induce ultrafast phase transitions into metastable collective states [START_REF] Li | Terahertz fieldinduced ferroelectricity in quantum paraelectric SrTiO 3[END_REF], for THz electron acceleration [START_REF] Nanni | Terahertz-driven linear electron acceleration[END_REF], ultrafast chemistry [START_REF] Larue | Thz-pulse-induced selective catalytic co oxidation on ru[END_REF] or even for atom probe tomography [START_REF] Vella | High-resolution terahertz-driven atom probe tomography[END_REF].

Terahertz time-domain spectroscopy [START_REF] Jepsen | Terahertz spectroscopy and imaging -modern techniques and applications[END_REF][START_REF] Xie | The application of terahertz spectroscopy to protein detection: A review[END_REF][START_REF] Neu | Tutorial: An introduction to terahertz time domain spectroscopy (thz-tds)[END_REF] (THz-TDS) is one of the most spread applications for pulsed THz radiation. THz-TDS is also used extensively from material science [START_REF] Hangyo | Terahertz time-domain spectroscopy of solids: A review[END_REF] to chemistry [START_REF] Fischer | Chemical recognition in terahertz time-domain spectroscopy and imaging[END_REF][START_REF] Bawuah | Advances in terahertz time-domain spectroscopy of pharmaceutical solids: A review[END_REF], so that several commercial systems of this important tool are already available on the market (see [START_REF] Neu | Tutorial: An introduction to terahertz time domain spectroscopy (thz-tds)[END_REF] and references therein). A major limitation of common THz-TDS systems is the bandwidth of the THz emitters and detectors, as well as the amplitude of the generated THz fields. The most common emitters and detectors are based on photoconductive antennas or optical rectification in nonlinear crystals. However, their bandwidths are usually below 10 THz due to absorption by optical phonons. Moreover, these solid-state devices are limited by the crystal sizes and their fragility to laser-induced damage.

Consequently, the further development of THz-TDS systems goes hand in hand with the search for suitable and robust THz emitters [START_REF] Zhang | Intense terahertz radiation: generation and application[END_REF][START_REF] Hafez | Intense terahertz radiation and their applications[END_REF]. An interesting alternative source is based on an air plasma generated by a two-color laser pulse [START_REF] Cook | Intense terahertz pulses by four-wave rectification in air[END_REF][START_REF] Kress | Terahertz-pulse generation by photoionization of air with laser pulses composed of both fundamental and secondharmonic waves[END_REF][START_REF] Kim | Terahertz emission from ultrafast ionizing air in symmetry-broken laser fields[END_REF]. The THz spectra emitted from such plasma are extremely broadband, with ranges up to 75 THz [START_REF] Kim | Coherent control of terahertz supercontinuum generation in ultrafast laser-gas interactions[END_REF] and even beyond [START_REF] Matsubara | Ultrabroadband coherent electric field from far infrared to 200 thz using air plasma induced by 10 fs pulses[END_REF]. Such a THz source is naturally ultrashort and can be straight-forwardly coupled with air biased coherent detection [START_REF] Karpowicz | Coherent heterodyne time-domain spectrometry covering the entire "terahertz gap[END_REF][START_REF] Ho | Design and performance of reflective terahertz airbiased-coherent-detection for time-domain spectroscopy[END_REF] (ABCD). The use of air as generation and detection medium has several advantages. First, broadband THz fields with high amplitude can be generated. Second, there is no risk of damaging the source or detection medium. The absence of solid materials also means that the occurrence of gaps in the THz spectra due to optical phonons is avoided.

In this paper we report on recent progress in the development and exploitation of a THz platform based on the above mentioned air plasma source and ABCD setup [START_REF] Prost | Air-photonics terahertz platform with versatile micro-controller based interface and data acquisition[END_REF]. In Sec. 2 the platform is briefly reviewed, with special emphasis on recent improvements of the signal-to-noise ratio. Rigorous numerical simulations of the whole platform, comprising the key optical elements shaping the input pump pulse, are reported in Sec. 3. Attention is paid to the numerical modeling of the ABCD scheme, from which we reveal and explain important differences between generated and measured THz spectra. By comparing to simplified pump configurations, we identify key parameters of the pump pulse determining the THz yield and spectral shape. Finally, in Sec. 4 we present the performances of the platform by measuring the absorption of D-phenylalanine and thymine samples. The resulting absorption spectra are successfully confronted with Fourier transform infrared (FTIR) measurements.

Optical setup

To produce a broadband THz radiation we use an air plasma created by a two-color ultra-short laser pulse [START_REF] Cook | Intense terahertz pulses by four-wave rectification in air[END_REF][START_REF] Kress | Terahertz-pulse generation by photoionization of air with laser pulses composed of both fundamental and secondharmonic waves[END_REF][START_REF] Kim | Terahertz emission from ultrafast ionizing air in symmetry-broken laser fields[END_REF]. The two-color pulse composed of fundamental (FH) and second harmonic (SH) ionizes the gaseous medium, here mostly dioxygen molecules. Acceleration of the free electrons by the two-color laser electric field leads to an asymmetric current containing a low frequency component that acts as radiation source in the THz spectral range [START_REF] Kim | Coherent control of terahertz supercontinuum generation in ultrafast laser-gas interactions[END_REF][START_REF] Babushkin | Ultrafast spatiotemporal dynamics of terahertz generation by ionizing two-color femtosecond pulses in gases[END_REF]. The detection of the THz radiation is performed by air biased coherent detection (ABCD) [START_REF] Karpowicz | Coherent heterodyne time-domain spectrometry covering the entire "terahertz gap[END_REF][START_REF] Dai | Detection of broadband terahertz waves with a laserinduced plasma in gases[END_REF][START_REF] Lu | Broadband terahertz detection with selected gases[END_REF]. This technique relies on four wave mixing in air between the THz field, a DC bias field and a FH probe field. The produced SH radiation is detected and contains the information about the THz field.

The optical setup presented in Fig. 1 is an improved version of the one used in Ref. [START_REF] Prost | Air-photonics terahertz platform with versatile micro-controller based interface and data acquisition[END_REF]. We use a Ti:Sapphire chirp-pulse amplifier that produces pulses of 30 fs at a repetition rate of 5 kHz with a central wavelength of 800 nm. Only 20 % of the total energy output of the amplifier is used for our experiments, which corresponds to 350 µJ. This initial FH beam is split between two arms, the first one (74%) is used for the THz generation and the second one (26%) for the ABCD detection. Compared to our previous work [START_REF] Prost | Air-photonics terahertz platform with versatile micro-controller based interface and data acquisition[END_REF] the energy ratio between the generation and detection arm has been adjusted. On one hand, the increase in energy on the first arm leads to an enhancement of the THz yield. On the other hand, lower energy for the detection prevents any plasma formation in the ABCD detection, thus making it more stable. To get rid of water vapor, which presents strong absorption lines in the THz range, the whole beam-path of the THz field is placed in a dry air enclosure (see gray area on Fig. 1).

To produce our two color laser pulse, we use second harmonic generation (SHG) in a 100 µm thick β-barium borate (BBO) crystal. This crystal is cut such that type I phase matching is optimized. At the output of the crystal, the SH pulse is orthogonally polarized with respect to the FH polarization. We use an ultra-thin half-waveplate at 800 nm (quartz, 45µm thickness) to turn the polarization axis of the FH pulse parallel to the SH polarization axis. This waveplate acts as an approximate full waveplate for the SH component, thus rendering the polarization of the SH pulse almost unchanged and only slightly elliptical. Both optics are placed before the geometrical focus of the lens, with the distance between BBO and focus being about 7 cm.

This two color laser pulse produces a plasma that is located a few millimeters before the geometrical focus of the lens. Oscillation of the free electrons inside the plasma gives rise to THz radiation that is collected and collimated by a gold coated off-axis parabolic mirror (OAP). A silicon wafer, placed in the collimated beam path, is used to filter out the two color optical pulse while transmitting 50% of the THz power [START_REF] Kaltenecker | Attenuation of thz beams: A "how to" tutorial[END_REF]. A second OAP is used to focus the THz field onto the sample. A second pair of OAPs collects and focuses the THz field for the detection between the high voltage electrodes of the ABCD setup.

The part of the laser pulse used as a probe for the ABCD scheme goes through a delay line, which allows to control its delay τ with respect to the incoming THz pulse on the fs time-scale. The probe pulse is then focused by a lens through a hole in the last OAP, such that both foci of THz and probe pulses overlap. Note that the half waveplate located in the pump arm, before the focusing lens, renders the polarization of THz pulse and probe pulse parallel. The high voltage electrodes are used to apply a static bias field in the same direction as the pulses polarization.

Our detection scheme relies on four wave mixing mediated by the third order nonlinearity in air. The mixing process produces a SH pulse and involves the FH probe field E ω and a low frequency field E Ω = E THz + E bias , composed of the sum of THz and bias field. We use a bandpass filter to isolate the SH pulse and detect it with an avalanche photodiode. The resulting SH signal S 2ω on the photodiode depends on the delay τ , and by neglecting the finite duration of the probe pulse, the signal can be expressed as [START_REF] Karpowicz | Coherent heterodyne time-domain spectrometry covering the entire "terahertz gap[END_REF] S 2ω (τ ) ∝ I 2 ω E 2 THz (τ ) + E 2 bias + 2E THz (τ )E bias , (1) with I ω being the intensity of the FH probe pulse.

To extract the THz field E THz (τ ) from the SH signal we perform two successive measurements, S + 2ω and S - 2ω , while applying static bias fields of opposite signs. The difference of those two measurements

S THz (τ ) = S + 2ω (τ ) -S - 2ω (τ ) ∝ I 2 ω E THz (τ )|E bias | (2
) is directly proportional to the THz field.

In this regards, the acquisition electronics has been improved compared to the one previously reported in [START_REF] Prost | Air-photonics terahertz platform with versatile micro-controller based interface and data acquisition[END_REF]. We adapted the time constant of our high voltage supplies, allowing us to increase the high voltage sign switch frequency while keeping the same stability on the high voltage output. We also replaced the 10 bits analog to digital converter integrated in our micro-controller by an external 15 bits one (ADS1115). This enhances the resolution and sensitivity of our detection significantly and gives us the opportunity to perform statistics on several measurements, while previously averaging over several measurements were performed to oversample the signal. This upgrade of the acquisition electronics coupled with the adjustment of the energy ratio between the pump and probe arms of the experiment gives a better overall stability with an improved resolution and sensitivity of the acquisition.

The measured THz pulses are almost single cycle with an envelope of about 50 fs FWHM duration. The measurements are performed with a bias field of E bias = ±10 kV/cm with steps of 5 fs over a span of 10 ps. By comparing the signal of the photodiode with and without bias field applied [START_REF] Iwaszczuk | Terahertz field enhancement to the mv/cm regime in a tapered parallel plate waveguide[END_REF], we measure a THz electric field strength of 1 MV/cm. Using the protocol described in [START_REF] Naftaly | Methodologies for determining the dynamic ranges and signalto-noise ratios of terahertz time-domain spectrometers[END_REF] we estimate the signal to noise ratio (SNR) and dynamic range (DR) of our acquisition to 350 and 10 4 , respectively. These estimations were performed in the temporal domain directly on the signals measured by the acquisition electronics, that is, S 2ω given by Eq. ( 1). Switching to a 15 bits analog to digital converter allows us to perform statistics on our acquisitions. By measuring the mean amplitude of the signal and its standard deviation, we calculate the SNR as the ratio between the two. The DR is calculated as the ratio between the maximum amplitude of the signal and the standard deviation of the noise signal. The noise signal is taken as the signal in absence of the THz field, in our case we measure it at long negative delays, that is, before the THz pulse arrival.

Numerical simulations

In this section, we present comprehensive numerical simulations of the different stages of our THz platform. Special efforts were made to meet the same field configurations as in the experiments, that is, propagation effects in all optical components after the focusing lens are taken into account. In particular, walk-off effects in the BBO crystal, in the dual wavelength plate as well as in the free air propagation turn out to be crucial for a proper characterization of the system. The air-plasma based THz source is modelled by a vectorial version of the unidirectional pulse propagation equation (UPPE) [START_REF] Kolesik | Unidirectional optical pulse propagation equation[END_REF][START_REF] Kolesik | Nonlinear optical pulse propagation simulation: From Maxwell's to unidirectional equations[END_REF] that governs the forward-propagating transverse electric field E = (E x , E y ) of elliptically-polarized pulses

∂ z ˆ E = i k 2 (ω) -k 2 x -k 2 y ˆ E + i µ 0 ω 2 2k(ω) ˆ P Kerr - µ 0 ω 2k(ω) ˆ J + ˆ J loss . (3)
We use such vectorial model to take into account the slightly elliptical polarization of the two-color pump induced by the ultra-thin waveplate. In Eq. (3), ˆ E(k x , k y , z, ω) is the Fourier transform of the transverse laser electric field with respect to x, y, and t. The longitudinal laser electric field component E z is neglected for the moderate focusing conditions used in the experiments. The square-root on the right-hand side of Eq. ( 3) describes linear dispersion and diffraction of the pulse, where k(ω) = ωn(ω)/c contains the frequency dependent refractive index n(ω) of air. The other terms contain the third-order nonlinear polarization

P Kerr (t) = 0 χ (3) E 2 (t) E(t) , (4)
with Kerr index n 2 = 3χ (3) /4n 2 0 c 0 . Dispersion of P Kerr is neglected, that is, χ (3) is taken as constant in the relevant frequency window, and consistently n 0 = n(ω FH ) is used in the prefactor. The electron current J follows from a Drude model

∂ t J(t) + ν c J(t) = q e m e N e (t) E(t) , (5)
and a loss term J loss (E, N e ) due to ionization [START_REF] Bergé | Optical ultrashort filaments in weakly-ionized, optically-transparent media[END_REF][START_REF] Bergé | Few-cycle light bullets created by femtosecond filaments[END_REF][START_REF] Babushkin | Ultrafast spatiotemporal dynamics of terahertz generation by ionizing two-color femtosecond pulses in gases[END_REF]. The free electron density N e is governed by

∂ t N e (t) = W (E)[N a -N e (t)] (6)
with field dependent ionization rate W (E) and density of neutral species N a to be ionized1 . In Eq. ( 5), q e and m e are the electron charge and mass, respectively, and ν c denotes the electronneutral collision rate. More details about this model can be found in [START_REF] Tailliez | Terahertz pulse generation by twocolor laser fields with circular polarization[END_REF] while the physical parameters for air are taken from [START_REF] Nguyen | Spectral dynamics of thz pulses generated by two-color laser filaments in air: the role of kerr nonlinearities and pump wavelength[END_REF].

The simulations employ the same laser parameters as our experiments, namely an incident linearly polarized 240 µJ laser pulse at 800 nm with 30 fs FWHM duration and 6 mm FWHM spatial beam width on the focusing lens L with f = 30 cm. This fundamental harmonic (FH) pulse is then frequency doubled in the 100 µm-thick BBO crystal situated 7 cm before the linear focus point. An orthogonally polarized co-propagating second harmonic (SH) pulse is created (5.5% energy conversion), where, in contrast to previous works [START_REF] Nguyen | Wavelength scaling of terahertz pulse energies delivered by twocolor air plasmas[END_REF], both spatial and temporal walk-off are taken into account by solving

i [∂ z + δ τ ∂ τ -δ x ∂ x ] A SH + χA 2 FH e i∆kz = 0 , (7) i∂ z A FH + χA * FH A SH e -i∆kz = 0 , (8)
for the slowly varying envelopes A SH and A FH expressed in a reference frame co-moving with the FH pulse (retarded time τ ). The symbol χ comprises the effective χ (2) nonlinearity, ∆k = 2k FHk SH is the phase mismatch, and

δ τ = 1/v SH g -1/v FH g
is the group velocity mismatch responsible for temporal walk-off. For type-I SHG the SH pulse is extraordinarily polarized, that is, both k SH and v SH g depend on the angle between propagation direction (wave vector) and optical axis of the uniaxial BBO crystal. This angle is chosen to 0.510 rad, such that the phase mismatch ∆k = 0 vanishes. Spatial walk-off is characterized by the walk-off angle δ x , that is, the angle between propagation direction (wave vector) and Poynting vector of the SH beam. Note that we defined the coordinate system such that spatial walk-off goes in x direction. For evaluating the relevant parameters δ τ and δ x we used the BBO dispersion relation given in Ref. [START_REF] Tamošauskas | Transmittance and phase matching of bbo crystal in the 3-5 µm range and its application for the characterization of mid-infrared laser pulses[END_REF], yielding δ τ = 1.93 × 10 -10 s/m and δ x = -0.678 rad.

After passing through the BBO crystal, the polarization of the FH and SH pulses is aligned by means of a 45 µm thick dual wavelength waveplate made of quartz. Because of the waveplate configuration, no spatial walk-off takes place. We ignore any nonlinear effects in the waveplate but take into account the temporal walk-off δ τ = 2.17 × 10 -10 s/m using the quartz dispersion relation given in Ref. [START_REF] Ghosh | Dispersion-equation coefficients for the refractive index and birefringence of calcite and quartz crystals[END_REF].

At a distance of about 4 mm before the geometrical focus, here denoted by z = 0, the twocolor pump intensity exceeds the threshold for ionization and an air-plasma is produced. Therefore, we start the computationally expensive UPPE simulation solving Eq. ( 3) at z = -5 mm. All other pump propagation after the lense, including SHG in the BBO, polarization rotation in the waveplate, and linear propagation in air is computationally light and handled by a pre-processor. We neglect the effect of the Kerr nonlinearity on the pulse propagation in the pre-processor. The value of the B-integral at z = -5 mm is smaller than one, and we checked that starting the UPPE simulation at z = -10 mm, which reduces the Bintegral by about a factor two, yields no significant differences. Figure 2(a) shows the resulting simulated peak electron density versus the propagation distance z. In the same figure, the energy buildup of the generated spectral components below 90 THz is plotted. Note that this electron density response and THz energy yield differ from those reported in [START_REF] Prost | Air-photonics terahertz platform with versatile micro-controller based interface and data acquisition[END_REF] because of a higher input beam energy and fraction of second harmonic. Interestingly, THz generation stops well before the end of the plasma zone. Before explaining this observation, it is worth noticing that the generated THz spectrum broadens significantly during the generation process, as shown in Fig. 2(b). Higher frequencies generated later along the propagation path, that is, farther away from the beginning of the plasma.

As already indicated above, the pre-processor takes into account the dispersion properties of the BBO, of the quartz dual wavelength waveplate as well as of the ambient air. This is crucial in order to get the correct beam and pulse profile as well as its (slightly elliptic) polarization state. In particular the temporal walk-off between the two pulse components turns out to be of utmost importance, in agreement with [START_REF] Zhang | Optimum chirp for efficient terahertz generation from two-color femtosecond pulses in air[END_REF]. To get further insight, the main FH and SH field components at the beginning (z = -3 mm) and at the end (z = -1 mm) of the region where THz radiation is generated are shown in Fig. 2(c)-(f). The orthogonal field components (not shown) are about one order of magnitude smaller. Obviously, the SH pulse is delayed by about 25 fs compared to the FH pulse, thus, the delay is of the order of the pulse duration itself. At z = -3 mm the two components overlap near zero delay, and the photocurrent mechanism for THz generation is active. At z = -1.5 mm, however, plasma defocusing in the trailing part of the FH pulse has completely separated the two colors in the spatio-temporal domain. This is the reason why THz generation stops roughly in the middle of the plasma at z = -1.5 mm, and the simulated THz energy yield remains of the order of few tens of nJ.

Figure 3(a) shows the spatial evolution of the THz fluency in the source area. During the main THz production stage z < -1.5 mm the THz beam is bell-shaped, but then it quickly evolves into the typical ring structure of conical emission [START_REF] Klarskov | Experimental threedimensional beam profiling and modeling of a terahertz beam generated from a two-color air plasma[END_REF]. The final THz spectral energy density at the end of the THz production stage (z = -1 mm) in Fig. 3(b) suggests that this broadband source could be used for spectroscopy in the range of 1 to 50 THz and beyond.

Finally, we want to shed some light on the THz detection process, that is, the part of the generated THz field that we can detect in our ABCD measurements. One could naively expect to measure a pattern close to the final spectrum at the end of the plasma in Fig. 3(b). This is, however, not the case. Both high and low THz frequencies are strongly suppressed in the ABCD scheme for the following two reasons.

First, the finite duration of the probe pulse (here 30 fs FWHM duration) poses an upper bound to the detectable THz frequencies [START_REF] Oh | Intense terahertz generation in two-color laser filamentation: energy scaling with terawatt laser systems[END_REF][START_REF] Andreeva | Ultrabroad terahertz spectrum generation from an air-based filament plasma[END_REF]. In the equation for the ABCD THz signal Eq. ( 2) the convolution with the temporal envelope of the probe pulse is omitted for simplicity. Taking this into account, one finds for the THz signal

S THz (τ ) ∝ |E bias | I 2 ω (t -τ )E THz (t)dt . ( 9
)
In temporal Fourier domain this implies that the THz spectrum ÊTHz (ω) gets multiplied by the Fourier transform of I 2 ω (t), that is, the longer the probe pulse duration, the more narrow the accessible spectral window. For our configuration, frequencies above 30 THz are strongly suppressed by the ABCD temporal optical resolution, as can be appreciated when comparing Figs. 4(b) and 3(b). Any features which affect this spectral region, e.g., the influence of a chirp in the pump pulse as reported in Ref. [START_REF] Prost | Air-photonics terahertz platform with versatile micro-controller based interface and data acquisition[END_REF], are thus not detectable in the present experimental configuration.

Second, the spectrum shown in Fig. 3(b) is spatially integrated over the transverse plane. In contrast to that, the ABCD setup detects only the part of the THz field overlapping with the focused IR probe pulse and the high voltage field. This is also reflected by Eq. ( 9) when one remembers that both I ω and E THz depend on spatial variables as well. Therefore, it is instructive to compute the THz field as it would appear in the detection area between the two electrodes. The four parabolic mirrors are essentially forming two consecutive 4f systems. Thus, an image of the THz source is produced. The resulting evolution of the THz fluence versus the propagation distance is shown in Fig. 4(a). Here, the 50% transmittance of the silicon wafer is taken into account. For convenience, the position of the image of the geometrical focus of the laser pump is set to z = 0, which makes it easier to compare with Fig. 3. The highest THz fluency (and intensity) in the detection area, as shown in Fig. 4(a), is obtained at z = -2.5 mm. This corresponds to the position of maximum THz generation in the source area, cf. Figs. 2 and3. As already mentioned, in order to perform ABCD measurements it is crucial to have a good overlap between THz and IR probe pulse. In the experiments, the probe pulse focus position is adjusted such that the ABCD signal is maximized. It is therefore reasonable to assume that the IR probe pulse is focused at the THz fluency maximum at z = -2.5 mm. The probe pulse focal FWHM spot size can be estimated to be less than 20 µm, that is, much smaller than the THz spot-size. Thus, only a small fraction of the THz field, close to the optical axis, actually contributes to the detection process. To account for this spatial overlap we modify the integral in Eq. ( 9) to

S THz (τ ) ∝ I 2 ω (t -τ, x, y)E THz (t, x, y)dtdxdy ,
where we assume that the THz signal S THz mainly originates from near the probe focal plane and thus ignore integration along z. Because the THz spot size is strongly frequency dependent, this leads to suppression of lower THz frequencies in the ABCD scheme. A similar effect was reported for THz detection via electro-optical sampling in [START_REF] Shipilo | Low-frequency content of thz emission from two-color femtosecond filament[END_REF]. Besides, parts of the low frequency field may not be collected by the parabolic mirror because beam divergence increases with decreasing frequency. This suppression of THz frequencies below 10 THz can be clearly seen in the simulated ABCD signal shown in Fig. 4(b), in particular when it is compared to the spectral energy density emitted by the source in Fig. 3(b).

The simulated ABCD spectrum is in excellent agreement with the experimental spectrum shown as a red dotted line in Fig. 4(b), obtained by taking the Fourier transform of the measured time dependent THz signal. The experimental spectrum exhibits a smooth profile extending from 1 THz to 30 THz, slightly narrower than the simulated one. We attribute this discrepancy to uncertainties in the laser pulse (probe) duration.

Even at this humidity level, water absorption lines can still be seen between 5 THz and 10 THz [49]. For this measurement, we used an average over 5 laser shots for each data point, which results in an acquisition time of about 10 minutes.

Before presenting THz-TDS results obtained with the platform it is worth to shed light on the impact of the pump configuration on both THz spectra and measured ABCD signals. To this end, we performed several simulations employing simplified synthetic two-color pump pulses as initial datum for the UPPE simulations starting at z = -5 mm as already mentioned above. The synthetic pump pulses neglect any previous propagation effects in the optical setup (except linear focusing), and the SH pulse profile is just proportional to the square of the FH pulse profile. We only allow two control parameters which turned out to have crucial impact on the THz generation, namely, the temporal separation between SH and FH pulse τ and phase-angle φ between SH and FH carrier wave. From the original pump field computed by the aforementioned pre-processor we extracted τ ≈ 25 fs and φ ≈ π/2 as closest match for the synthetic pump field at z = -5 mm (see also pump snapshots shown in Fig. 2). The resulting THz spectrum in the source area is shown as the black dashed line in Fig. 5(a). Obviously, the THz spectral shape is not very different compared to the original one (solid blue line), but the yield is reduced by roughly a factor five. We checked that this difference in the THz yield is due to a slightly larger SH pulse duration and thus a better overlap for the original pump configuration. The similarity of the spectral shape however comes from the matching of the phase angle φ. The same synthetic pulse but with a phase angle φ = 0 produces a quite distinct spectrum, as can be inferred from the black dotted line in Fig. 5(a). By contrast, removing the temporal separation by setting τ = 0 increases the THz yield by more than one order of magnitude and enlarges the spectral extend (black dashed-dotted line), underlining once more the importance of this parameter. Other features of the original pump profile, like deviation from linear polarization and slight spatial misalignment can be seen in the THz field profiles (not shown) but are of minor influence on the THz yield and spectral energy density.

Finally, it is worth checking if the above modifications of the THz spectra would be detectable by the ABCD scheme on our platform. Figure 5(b) shows the resulting simulated signals for our 30 fs probe pulse. Color codings are the same as in Fig. 5(a). As we have already seen above, the high frequency part > 30 THz is completely lost, so that the spectral shapes of all curves become very similar. Some variations due to changes in the phase angle survive, however, they are not very significant. Merely the significant differences in the THz yield are clearly visible in the ABCD THz signal, as expected. Using significantly shorter probe pulses than the 30 fs employed in our setup would render the ABCD method capable of detecting those shape differences.

Time-domain spectroscopy

The THz platform allows to perform THz-TDS by placing a sample at the intermediate focus in the THz beam path. To demonstrate this capability of the improved setup, we use D-phenylalanine and thymine samples. Crystalline powders were obtained from Sigma-Aldrich, with purities of 98% and 99%, respectively. Our samples have been prepared by deposition of the molecules on small pieces of silicon wafer (∼1 cm 2 ). This deposition was performed by sublimating some powder under vacuum, while maintaining the temperature below the melting point.

For those measurements, the temporal stepsize in the ABCD is reduced to 2 fs to improve the resolution of the signal in temporal domain, thus expanding the spectral span to 250 THz, and the averaging is performed over 10 laser shots per data point. To keep the acquisition time below 30 minutes the temporal span is reduced to 5 ps, which increases the spectral step-size to 0.2 THz. Figure 6(a) shows the resulting temporal profile of the THz electric field for D-phenylalanine, together with a reference measurement of a clean piece of silicon wafer. The strong modification of the THz pulse by the propagation through the deposited D-phenylalanine is clearly visible. The most noticeable change is the (small) delay of the main pulse caused by the different group velocity inside the sample. Furthermore, a stretching of the signal induced by absorption and dispersion in the sample leads to an oscillation like pattern, as can be seen between 0.1 and 0.6 ps. Once transformed to the spectral domain, those temporal modifications result in absorption bands between 5 THz and 30 THz, see Fig. 6(b). As in the reference measurement shown in Fig. 4(b), the relative humidity was kept below 5%, nevertheless, spectral lines attributed to water can still be observed below 10 THz. The absorption line at 18 THz appears due to a phonon overtone in the silicon wafer [START_REF] Kaltenecker | Attenuation of thz beams: A "how to" tutorial[END_REF].

To distinguish the absorption lines induced by the molecular samples from the ones coming from environment and substrate, we normalize the transmission spectra by dividing the measured

Conclusion

In this paper, we have presented a detailed characterization of an air photonics based terahertz platform and its ability to perform absorption spectroscopy in bio-relevant molecular samples. The whole platform, including the detection scheme, was carefully benchmarked against comprehensive numerical simulations. This confrontation of theoretical results and measurements performed on our setup revealed that it is important to take into account pulse and beam shaping effects in the setup even before the plasma for THz generation is formed. In particular, walk-off effects in the SHG process as well as along the full beam path (e.g., waveplate and two-color propagation in air) turned out to be crucial, because they alter the spatio-temporal overlap of the two pump colors in the plasma significantly. Moreover, we showed that the THz spectra retrieved by ABCD measurements are strongly dependent on the probe characteristics, namely its spatial and temporal extends. While the latter determines the high-frequency cut-off of the measurable THz spectrum, the former one sets a limit on the low frequency part.

With the upgraded precision of the acquisition electronics, we report a signal to noise ratio of above 300 and a dynamic range of 10 4 , in the temporal domain, which makes our platform suitable for THz time-domain spectroscopy. We demonstrate this capability by measuring transmission spectra of D-phenylalanine and Thymine samples. Because of its simple layout, other types of applications employing THz pulses as pump or probe such as THz streaking [START_REF] Ardana-Lamas | Temporal characterization of individual harmonics of an attosecond pulse train by thz streaking[END_REF], molecular orientation [START_REF] Xu | Longlasting molecular orientation induced by a single terahertz pulse[END_REF] or nonlinear spectroscopy [START_REF] Tcypkin | Giant third-order nonlinear response of liquids at terahertz frequencies[END_REF] are possible.

Fig. 1

 1 Fig. 1 Schematic of the THz platform. BS : Beam splitter (76% reflection); L : focusing lens (f = 300 mm); DL : Delay line; HWP : Half waveplate; BH : composed of a BBO crystal and a dual wavelength waveplate; BD : Beam dump; SiW : Silicon wafer; S : Sample; HV : High voltage electrodes; F : 400 nm band-pass filter; APD : Avalanche photodiode. The area in gray indicates the dry air enclosure.

  Spectral Energy Density [arb.u.] (b) z = -1.0 mm z = -2.0 mm z = -2.5 mm z = -3.0 mm z = -3.5 mm z = -4.0 mm FH @ z = -3 mm, y = 0 µm

Fig. 2

 2 Fig. 2 THz source: (a) Peak plasma density and THz energy (ν < 90 THz) versus propagation distance relative to the geometrical focus at z = 0. (b) Snapshots of the transverse integrated THz spectrum at several positions z along the plasma. Snapshots of the main FH (c) and SH (d) polarization component of the pump field near the beginning of the plasma at z = -3 mm. (e) and (f) show the same fields inside the plasma after THz generation has stopped due to lack of temporal overlap between SH and FH component at z = -1 mm.

Fig. 3

 3 Fig. 3 THz source: (a) THz fluency (ν < 90 THz) in the source area as a function of the propagation distance. The position of the geometrical focus of the pump laser is denoted with z = 0. (b) Transverse integrated THz spectrum at z = -1 mm, that is, at the end of the zone of THz generation [same as the blue curve in Fig. 2(b), but in linear scale].

Fig. 4

 4 Fig. 4 THz detection: (a) THz fluency (ν < 90 THz) in the detection area between the HV electrodes as a function of the propagation distance. The position of the image of the geometrical focus of the pump laser is denoted with z = 0. (b) Spectrum from a simulated ABCD method involving a 30 fs probe pulse focused to 15 µm waist at z = -2.5 mm (see text). The experimental spectrum taken at a relative humidity below 5% inside the enclosure is shown as red dotted line.

1 Fig. 5

 15 Fig. 5 Impact of pump distortions on the THz spectra in (a) the source area and (b) on the ABCD signal. The blue lines show the spectra as reported in Fig. 3(b) and Fig. 4(b) but in logarithmic scale, respectively. The black lines show what happens when the pump pulse at z = -5 mm is replaced by a synthetic two-color pulse with delay τ and phase-angle φ as indicated in the legend (see text for details).

Fig. 6

 6 Fig. 6 (a) Temporal profile of the THz electric field with a sample of D-phenylalanine deposited on a silicon wafer (blue curve) and without deposition on the silicon wafer as reference (red curve). Both spectra are normalized to the reference spectrum peak value. (b) Retrieved intensity spectra with the same color code.

Fig. 7

 7 Fig. 7 Reconstructed transmission spectrum with (a) Dphenylalanine and (b) thymine samples. The blue lines correspond to the spectra measured by THz-TDS and the red ones to measurements done with a FTIR spectrometer. The vertical dashed lines are a guide for the eye indicating the main line positions.

For our laser configuration only oxygen molecules, which have the lowest ionization potential of main air constituents, get ionized.

Acknowledgments. We thank B. Moge, I. Aguili, F. Khalid Balyos and C. Clavier for technical support. We thank P. U. Jepsen and B. Zhou from the Technical University of Danemark for fruitful discussions. We acknowledge financial support from CNRS and the ANR ALTESSE2 (ANR-19-ASMA-0007) project. S. Skupin thanks the Qatar National Research Fund (NPRP 12S-0205-190047) for support. Numerical simulations were performed using resources at Grand Équipement National De Calcul Intensif (GENCI) (A0100507594).

Statements and Declarations

Data Availability Statement

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Conflict of Interest

The authors have no conflicts to disclose.