Mathias Rapacioli
email: mathias.rapacioli@irsamc.ups-tlse.fr

Nathalie Tarrat
email: nathalie.tarrat@cemes.fr

Periodic DFTB for supported clusters : Implementation and application to benzene dimers deposited on graphene

Keywords: Periodic DFTB, Benzene dimers, Graphene 14

The interest for properties of clusters deposited on surfaces has grown in the recent 1 years. In this framework the Density Functional based Tight Binding (DFTB) method appears as a 2 promising tool due to its ability to treat extended systems at the quantum level with a low com-3 putational cost. We report the implementation of periodic boundary conditions for DFTB within 4 the deMonNano code with k-points formalism and corrections for intermolecular interactions. 5 The quality of the DFTB calculations is evaluated by comparison with dispersion-corrected DFT 6 calculations. Optimized lattice properties for a graphene sheet and graphite bulk are in agreement 7 with reference data. The deposition of both benzene monomer and dimers on graphene are inves-8 tigated and the observed trends are similar at the DFT and DFTB levels. Moreover, interaction 9 energies are of similar orders of magnitude for these two levels of calculation. This study has 10 evidenced the high stability of a structure made of two benzene molecules deposited close to each 11 other on the graphene sheet. This work demonstrates the ability of the new implementation to 12 investigate surface deposited molecular clusters properties.

13

Introduction 15

The modelling of functional extended surfaces has grown in past decades to in-16 vestigate, for fundamental and engineering purposes, a large number of phenomena 17 or applications such as e.g. deposition [START_REF] Sadowska | Nanoparticle deposition on heterogeneous sur-313 faces: Random sequential adsorption modeling and experiments[END_REF], growth and migration [START_REF] Li | Modeling the Migration of Platinum Nanoparticles on Surfaces Using a Kinetic 317 Monte Carlo Approach[END_REF], 2D assembly [START_REF] Rochefort | Role of the Structure and Reactivity of Cu and 321 Ag Surfaces in the Formation of a 2D Metal-Hexahydroxytriphenylene Network[END_REF],

18 catalysis [START_REF] Bruix | First-principles-based multiscale modelling of 324 heterogeneous catalysis[END_REF], electrocatalysis [START_REF] Abidi | Atomistic modeling of electrocatalysis: Are 326 we there yet[END_REF], photocatalysis [START_REF] Zhang | Surface modeling of photocatalytic materials for water 328 splitting[END_REF], molecular electronics [START_REF] Vilan | Chemical Modification of Semiconductor Surfaces for Molec-330 ular Electronics[END_REF], depollution 19 [START_REF] Shahmoradi | Removal of methylmercaptan pollution using Ni and Pt-decorated 334 graphene: an ab-initio DFT study[END_REF], sensing [START_REF] Li | DFT calculation and analysis of the gas sensing 337 mechanism of methoxy propanol on Ag decorated SnO2 (110) surface[END_REF]... Many of these studies have focused on deposited clusters, i.e. finite 20 aggregations of basis elements (atoms or molecules) adsorbed on surfaces. Indeed, the 21 physico-chemical properties of a cluster are distinct from the ones of both the single 22 entities and the infinite cluster (bulk) and strongly depend on the size and structure of 23 the cluster [START_REF] Soini | Size-dependent properties of transition metal clusters: from molecules 340 to crystals and surfaces -computational studies with the program ParaGauss[END_REF]. Understanding and controlling the structuration of deposited clusters 24 could thus allow a precise tuning of their properties.

25

The theoretical study of clusters deposited on extended surfaces is very challenging 26 due to the size of the space of structural and electronic configurations to be explored [START_REF] Elstner | Self Consistent-Charge Density-391 Functional Tight-Binding Method for Simulations of Biological Molecules[END_REF] and to the high level of computational methods that has to be implemented. Indeed, 28 the a priori unknown nature of the interactions between the cluster building blocks 29 and between the cluster and the surface (with potential charge transfers at the cluster- The Density Functional based Tight-Binding method (DFTB) can be derived from 67 DFT from several approximations [START_REF] Porezag | Construction of tight-343 binding-like potentials on the basis of density-functional theory: Application to carbon[END_REF][START_REF] Elstner | Density functional tight binding[END_REF][START_REF] Spiegelman | Density-functional tight-binding: basic concepts and applications to molecules and 353 clusters[END_REF][START_REF] Koskinen | Density-Functional Tight-Binding for Beginners[END_REF][START_REF] Elstner | Self Consistent-Charge Density-391 Functional Tight-Binding Method for Simulations of Biological Molecules[END_REF].

68

The first one relies on expression of molecular orbitals (MOs) φ i (r) as linear combinations of atomic orbital (LCAO)-type basis sets using minimal valence bases χ µ .

φ i (r) = ∑ µ c iµ χ µ (r) (1)
A Taylor expansion of the DFT energy is done as a function of the electronic density, the real density ρ of the system minimizing the Kohn-Sham energy being searched as a perturbation with respect to a reference density ρ 0 (ρ = ρ 0 + δρ) :

E[ρ(r)] = E[ρ 0 (r)]+ δE[ρ(r)] δρ(r) ρ 0 δρ(r) + 1 2 δ 2 E[ρ(r)] δρ(r)δρ(r ) ρ 0 δρ(r)δρ(r )+ ... + 1 p! ... δ p E[ρ(r)] δρ(r)δρ(r )...δρ(r (p) ) ρ 0 δρ(r)δρ(r )....δρ(r (p) ) (2)
In the original version, also known as the non self-consistent DFTB (sometimes 69 referred to as zeroth-order DFTB or simply DFTB [START_REF] Porezag | Construction of tight-343 binding-like potentials on the basis of density-functional theory: Application to carbon[END_REF]12]), only the zeroth and first order 70 terms of the Taylor expansion are retained. In the DFTB2 scheme [START_REF] Elstner | Self Consistent-Charge Density-391 Functional Tight-Binding Method for Simulations of Biological Molecules[END_REF], also known as 71 self-consistent charge (SCC) DFTB, and in the DFTB3 scheme [START_REF] Yang | Extension of the self-consistent-charge density-394 functional tight-binding method: third-order expansion of the density functional theory 395 total energy and introduction of a modified effective coulomb interaction[END_REF], the second and third 72 order terms are also taken into account, respectively.

At the DFTB0 level, the potential energy reads :

E DFTB0 = ∑ α<β E rep (r αβ ) + ∑ iµν n i c iµ c iν H 0 µν (3)
with E rep (r αβ ) a repulsive contribution between atoms α and β, n i the occupation of 74 the orbital i and H 0 µν the matrix elements associated to the Kohn Sham operator at the 75 reference density expressed in the atomic basis. Its matrix elements, as well as those 76 of the atomic overlap matrix S, can be parameterized as only one-or two-body terms.

77

This is allowed by the definition of the reference density as a superposition of atomic Focusing from now on the SCC-DFTB level [START_REF] Elstner | Self Consistent-Charge Density-391 Functional Tight-Binding Method for Simulations of Biological Molecules[END_REF], the previous energy expression becomes :

78 densities ρ 0 = ρ α 0 + ρ β 0 + ρ γ 0 + ...
E SCC-DFTB = E DFTB0 + 1 2 ∑ αβ γ αβ q α q β (4)
The last term corresponds to the second-order contribution and depends on the electronic density fluctuation δρ represented by atomic charges q α . γ αβ is a matrix whose diagonal terms are equal to the atomic Hubbard parameters and off-diagonal terms contain the 1/R coulomb interaction between atomic charges and an exchange-correlation energy contribution:

γ αβ = 1 |r -r | + δ 2 E xc δρ(r)δρ(r ) ρ 0 F α 0 (r -r α )F β 0 (r -r β )drdr (5)
where F α 0 is the normalised spatial extension for the excess/default of electrons around 84 atom α with respect to the neutral atom, assumed here to have no angular dependence.

85

Since the second-order term contains atomic charges, this introduces a term depending on the charges H 1 (q) into the TB operator.

(H 0 + H 1 (q))C i = i SC i (6)
with

H 1 µν = 1 2 S µν ∑ ξ q ξ (γ αξ + γ βξ )
where µ and ν belong to atoms α and β, respectively. As the charges depend on the MO 86 coefficients c iµ , the new secular equation must be solved self-consistently with respect to 87 atomic charges, at the origin of the method's name self-consistent-charge (SCC-)DFTB.

88

In the standard SCC-DFTB version [START_REF] Elstner | Self Consistent-Charge Density-391 Functional Tight-Binding Method for Simulations of Biological Molecules[END_REF], the atomic charges are computed from the density matrix P and the atomic basis overlap S matrix within the Mulliken approximation.

q α = ∑ µ∈α ∑ ν P µν S µν ( 7
)
with

P µν = ∑ µν n i c iµ c iν
In previous works, we have shown that atomic charges can be improved by taking into account the bond polarisation, adapting the Charge Model class IV scheme for DFTB [START_REF] Rapacioli | Correction for dispersion and Coulombic interactions in molecular clusters with density 388 functional derived methods: Application to polycyclic aromatic hydrocarbon clusters[END_REF][START_REF] Li | New Class IV Charge Model for Extracting Accurate 398 Partial Charges from Wave Functions[END_REF]. This approach, which requires the calculation of Mayer's bond order, is computationally expensive and hardly transferable to a periodic implementation. To circumvent this bottleneck, we have recently introduced a simpler scheme, named in the following WMull for Weighted Mulliken charges [START_REF] Michoulier | Theoretical determination of adsorption and ionisation energies of polycyclic aromatic 402 hydrocarbons on water ice[END_REF], to correct atomic charges with the following expression

q α = ∑ µ∈α ∑ ν P µν S µν (1 + t αβ ) (8)
t αβ = -t βα is an empirical parameter accounting for a non-symmetric repartition of the electrons between different atomic types, the Mulliken symmetric repartition being recovered for t αβ = 0. The second order contribution to the Kohn Sham operator matrix is modified as follows :

H 1 µν = 1 2 S µν ∑ ξ q ξ (γ αξ (1 + t αβ ) + γ βξ (1 -t αβ )) (9)
We have shown that this simple scheme provides similar results to those obtained with 89 the Charge Model approach to model clusters of PAHs [START_REF] Dontot | Dependence upon charge of the 405 vibrational spectra of small Polycyclic Aromatic Hydrocarbon clusters: the example of 406 pyrene[END_REF] and water [START_REF] Simon | Contribution 408 of the density-functional-based tight-binding scheme to the description of water clusters: 409 methods, applications and extension to bulk systems[END_REF]. The former implementation of periodic DFTB within deMonNano was restricted to 92 the Γ-point approximation only [22]. In the present implementation, the electronic prob- For a given k-point, the molecular orbitals φ k i must full-fill the Bloch-Theorem, that is TR

φ k i = e ikR φ k i (r) ( 10
)
where TR is the operator associated to a translation of R, where R is a vector connecting two unit cells. This is achieved by expanding the MOs on an basis of Bloch functions built from the real space atomic orbitals.

φ k i (r) = ∑ k c k iµ χ k µ (r) with χ k µ (r) = 1 √ N ∑ N e ikR N χ µ (r -R N ) (11)
where the infinite sum relies on all the possible translation from the main unit cell to any of the other ones. The overlap and Kohn Sham operator matrices expressed in this basis can be written from their real-space equivalent making use of the following transformation :

H k µν = ∑ N e ikR N H µν (r α -r β -R N ) = ∑ N e ikR N (H 0 µν (r α -r β -R N ) + H 1 µν (r α -r β -R N )) and S k µν = ∑ N e ikR N S µν (r α -r β -R N ) (12)
where µ and ν belong to atoms α and β, respectively. 100

In the previous expressions, the matrix elements of H 0 and S are easily obtained from the DFTB Slater Koster tables and rapidly vanish for large values of R N . The first order contribution to H µν is also short range with respect to (r αr β -R N ) but contains an infinite long range coulomb sum :

H 1,k µν (r α -r β -R N ) = 1 2 S µν (r α -r β -R N ) ∑ ξ ∑ N q ξ (γ αξ (r α -r ξ -R N ) + γ βξ (r β -r ξ -R N )) (13)
In practice, this infinite sum is replaced by an Ewald summation. The secular equation is solved for each k-point :

H k C k i = k i S k C k i ( 14
)
The eigenvalues k i resulting from all the k-point secular equations are then sorted in ascending order to drive the determination of n k i , the orbital occupation number following either a canonical occupation or a Fermi distribution.

The density matrix can therefore be computed for each k-point

P µν (k) = ∑ i n k i c k * iµ c k iν ( 15
)
We follow the approach of reference [START_REF] Aradi | a Sparse Matrix-Based Implementation of 355 the DFTB Method[END_REF], which consists in building the real space density matrix and computing atomic charges in the real space. The real space density matrix is obtained by summing over the k-points :

P µν (R N ) = ∑ k P µν (k)e -ikR N (16)
and Mulliken charges are then computed as follows :

q α = ∑ R N P µν (R N )S µν (r α -r β -R N ) (17)
and used as inputs for the next SCC cycle.

The simple WMull correction to Mulliken charges can be generalised to the periodic equations replacing equation 17 and 13 by

q α = ∑ R N P µν (R N )S µν (r α -r β -R N )(1 + t αβ ) (18)
and

H 1,k µν (r α -r β -R N ) = 1 2 S µν (r α -r β -R N ) ∑ ξ ∑ N q ξ (γ αξ (r α -r ξ -R N )(1 + t αβ ) + γ βξ (r β -r ξ -R N )(1 -t αβ ))
2. Computational details

DFTB calculations

Different DFTB parameters are available in the literature (cf. website www.dftb.org), depending on the choices made during the parameterization procedure such as the DFT functional, the basis sets type (Gaussian, Lorentzian) used to generate the atomic orbitals, the confinement imposed on these orbitals, the reference data used to compute the repulsive contribution E rep and, for the second and third order DFTB, the values of the atomic Hubbard parameters and their derivatives. In this work, we are working with the BIO DFTB set of parameters [START_REF] Elstner | Self Consistent-Charge Density-391 Functional Tight-Binding Method for Simulations of Biological Molecules[END_REF] provided within the deMonNano code (equivalent to the mio parameters from the website www.dftb.org). Dispersion interaction corrections can be introduced in the DFTB Hamiltonian using an empirical diatomic formulae. Two types of corrections are available in the deMonNano code and will be tested in the next section. The first one (hereafter labelled D1) is a Lennard Jones type potential with short range corrections introduced by Zhechkov et al. [START_REF] Zhechkov | An Efficient a Posteriori 412 Treatment for Dispersion Interaction in Density-Functional-Based Tight Binding[END_REF]. The second one (hereafter labelled D2 [START_REF] Rapacioli | Correction for dispersion and Coulombic interactions in molecular clusters with density 388 functional derived methods: Application to polycyclic aromatic hydrocarbon clusters[END_REF]) is given by

E disp = -∑ N ∑ α,β f damp (|r α -r β -R N |) C 6 AB |r α -r β -R N | 6 (19)
where f damp is a damping function screening the short range contribution and C 6 αβ 111 an empirical parameter (see [START_REF] Rapacioli | Correction for dispersion and Coulombic interactions in molecular clusters with density 388 functional derived methods: Application to polycyclic aromatic hydrocarbon clusters[END_REF] for details). In both cases, only the van der Waals 112 contributions larger that 10 PAW pseudopotentials [START_REF] Blöchl | Projector augmented-wave method[END_REF][START_REF] Kresse | From ultrasoft pseudopotentials to the projector augmented-wave 425 method[END_REF] and the DFT-D3 semiempirical dispersion-corrected func-124 tional in its zero-damping formalism [START_REF] Grimme | A consistent and accurate ab initio parametriza-427 tion of density functional dispersion correction (DFT-D) for the 94 elements H-Pu[END_REF]. This functional has been chosen as it has been 125 reported as a relevant choice for studies involving graphene [START_REF] Lebedeva | Comparison of performance 430 of van der Waals-corrected exchange-correlation functionals for interlayer interaction in 431 graphene and hexagonal boron nitride[END_REF]. A conjugate-gradient 126 algorithm was used to relax the ions and the convergence criterion was set up so that 127 the maximum atomic force was less than 0.01 eV Å -1 , all atoms being allowed to relax 128 unconstrained. To avoid interactions between the benzene monomers/dimers and their 129 periodic images, a cubic box measuring 50Å on a side was used for isolated systems. 1), with one k-point in 165 the x and y directions. The total energy varies by 2 10 -4 eV/atom (0.0004%) when the 166 number of k-points is increased from one to three in the z direction, and then remains 167 constant for calculations performed with five, seven and nine k-points in the z-direction.

168

We have thus determined the equilibrium parameters of the graphite bulk with three 169 k-points in the z directions (Table 1). Using either one or three k-points in the x and have been probed, labelled a 1 , a 2 , a 3 , a 4 , which can be visualised in Figure 1. The three 186 first ones correspond to the structures labelled a 1 , a 2 , a 3 in reference [START_REF] Zhang | Tuning the electronic 449 structure and transport properties of graphene by noncovalent functionalization: effects of 450 organic donor, acceptor and metal atoms[END_REF]; hollow, bridge, 187 top in reference [START_REF] Berland | Analysis of van der Waals density functional components: Binding 453 and corrugation of benzene and C 60 on boron nitride and graphene[END_REF] and AA, SP and AB in reference [START_REF] Ershova | Study of polycyclic aromatic hydrocarbons 456 adsorbed on graphene using density functional theory with empirical dispersion correction[END_REF]. The last structure a 4 was 188 named top-rot in reference [START_REF] Berland | Analysis of van der Waals density functional components: Binding 453 and corrugation of benzene and C 60 on boron nitride and graphene[END_REF] and also studied in reference [START_REF] Otyepková | Organic adsorbates have higher 459 affinities to fluorographene than to graphene[END_REF]. Table 2. Binding energies of benzene on graphene in eV. The experimental binding energy of a benzene molecule on a graphite surface is -0.50±0.08 eV [START_REF] Zacharia | Interlayer cohesive energy of graphite from ther-462 mal desorption of polyaromatic hydrocarbons[END_REF]. -0.428 -0.450 -0.453 -0.450 LDA [START_REF] Zhang | Tuning the electronic 449 structure and transport properties of graphene by noncovalent functionalization: effects of 450 organic donor, acceptor and metal atoms[END_REF] -0.16 -0.23 -0.24 ωB97X-D [START_REF] Ershova | Study of polycyclic aromatic hydrocarbons 456 adsorbed on graphene using density functional theory with empirical dispersion correction[END_REF] -0.47 optB86b-vdw [START_REF] Otyepková | Organic adsorbates have higher 459 affinities to fluorographene than to graphene[END_REF] -0.5 vdW-DF1 [START_REF] Berland | Analysis of van der Waals density functional components: Binding 453 and corrugation of benzene and C 60 on boron nitride and graphene[END_REF] -0.49 vdW-DF2 [START_REF] Berland | Analysis of van der Waals density functional components: Binding 453 and corrugation of benzene and C 60 on boron nitride and graphene[END_REF] -0.43 Expt. Saturated Adsorption Enthalpy [START_REF] Otyepková | Organic adsorbates have higher 459 affinities to fluorographene than to graphene[END_REF] -0.5 DFTB-D1 and DFTB-D2 results agree on the main trends, also present at the DFT-D3 [START_REF] Berland | Analysis of van der Waals density functional components: Binding 453 and corrugation of benzene and C 60 on boron nitride and graphene[END_REF] 3.6 vdW-DF2 [START_REF] Berland | Analysis of van der Waals density functional components: Binding 453 and corrugation of benzene and C 60 on boron nitride and graphene[END_REF] 3.5

DFTB-D1

Method

Benzene dimers in vacuum

201

Reproducing the benzene dimer potential energy surface is a challenging task for 202 DFT schemes, due to the fine competition between the various contributions to the total 203 energy. This is even more true for approximated schemes like the DFTB method. Briefly, The binding energies for the different optimized structures are reported in Table 4 217 for both the DFTB, DFT-D3 and ab initio reference calculations (CCSD(T) and SAPT). In Csoa -0.146 -0.152 -0.106 T -0.138 CCSD(T) [START_REF] Podeszwa | Potential Energy Surface for the Benzene Dimer 384 and Perturbational Analysis of pi-pi Interactions[END_REF] Csoa -0.12 -0.12 -0.07 SAPT [START_REF] Podeszwa | Potential Energy Surface for the Benzene Dimer 384 and Perturbational Analysis of pi-pi Interactions[END_REF] Csoa -0.12 -0.12 -0.08 shows that the latter is the most stable investigated configuration as it maximizes the 277 interaction between each benzene molecule and the graphene surface while maintaining 278 some stabilizing intermolecular interactions between the two benzene units. Again, this 279 conclusion holds at both the DFT and DFTB levels. It should also be noted that the 280 values of the interaction energies are of similar order for these two levels of calculation.

30

  surface interface) prevents the use of empirical force fields in favor of methods in which 31 the electronic structure is explicitly considered. The very high computational cost of 32 ab initio and Density Functional Theory (DFT) methods prohibiting their use for the 33 study of such systems, particularly when dealing with global optimisation or finite-34 temperature molecular dynamics, one strategy consists in implementing in a periodic 35 formulation approximate quantum mechanical methods. Among them, DFTB is an 36 approximated DFT schemes with a much lower computational cost enabled by the use of

90 1 . 2 .

 12 DFTB for periodic systems 91

93

  lem is searched self consistently after defining a set of k-points in the reciprocal space. A 94 step of the self-consistent scheme consists in solving separately the secular equations 95 for each k-point to obtain the molecular orbitals φ k i . Molecular orbitals obtained for all 96 k-point are then used to build the total electronic density, the latter being used as an 97 input for the next self-consistent step.98 99

130 3 .

 3 For supported ones, the calculations were performed on a 29.92Å x 34.55Å graphene 131 surface (these values having been calculated on the basis of the graphene equilibrium 132 lattice parameter reported in section 3.1) containing 392 carbon atoms placed in a 50 133 Å high simulation box to avoid any interaction between the adsorbed molecule and 134 the underside of the graphene sheet of the upper periodic box. Since the size of the 135 supercell was large enough, the Brillouin zone sampling in reciprocal space restricted 136 to the Γ-point was sufficient to ensure good convergence of the total energy, except 137 for the calculations aiming at determining the equilibrium parameters of the graphite 138 bulk which required a 1x1x5 k-points grid. A plane-wave kinetic energy cut-off of 450 139 eV was employed. For dealing with the partial occupancies around the Fermi level, a 140 Methfessel-Paxton smearing was used with σ = 0.2 eV [41]. 141 142 Results and discussion 143 In this section, we discuss the results of the DFTB calculations. All details about 144 the dispersion-corrected DFT calculation performed to evaluate the quality of the DFTB 145 calculations are given in the previous section, so that the computational details reported 146 below only concern the DFTB calculations.

147 3 . 1 .

 31 Graphene and Graphite148In order to model graphene, we have first optimized the lattice parameter, working 149 with a periodic box containing 392 atoms (∼30Å x 35Å x 50Å). For such a large simulation 150 box, the Γ-point approximation remains valid as the energy varies by less than 3.4 10 -4 151 eV/atom (0.0008%) when going from one to three k-points in the x and y directions, and 152 by less than 1.3 10 -4 eV/atom (0.0003%) when going from three to five k-points. The 153 equilibrium C-C bond length determined with one or three k-points in x and y directions 154 are the same at the precision of 10

170yTable 1 .

 1 directions led to the sames results at the target precision of 0.001 Å. The DFTB-D1 171 and DFTB-D2 C-C bond length are reduced by 0.001 Å with respect to their values in 172 the graphene sheet, a trend also observed at the DFT-D3 level. The DFTB-D1 interlayer 173 equilibrium distance (3.383±0.001 Å) is in agreement with both the value of reference 174 [33] with a four layers model in the Γ-point approximation (3.38 Å) and the experimental 175 data (3.356 Å). The interlayer distance is reduced to 3.131±0.001 Å at the DFTB-D2 level. 176 With respect to theoretical references (DFT-D3, RPA and QMC) and experimental values, 177 we can conclude that the DFTB-D1 method gives better quality results for the graphite 178 interlayer distances, while the the DFTB-D2 method prevails for the C-C bond length. 179 Graphene and Graphite structural data (in Å). *In these computational studies the d C-C distance was fixed to the one determined experimentally. * * Differences between the two DFTB-D1 calculations are detailed in the text.

Figure 1 .

 1 Figure 1. Isolated benzene molecule deposited on graphene.

  190 level : three almost degenerated structures, namely a 2 , a 3 and a 4 , and the a 1 structure 191 being less stable by about 0.012(DFTB-D2)/0.015(DFTB-D1)/0.025(DFT-D3) eV (see

204

  three characteristic structural families can be identified, namely sandwich (S), parallel-205 displaced (PD) and T-shaped (T), each one presenting several minima. In this work, the 206 sandwich eclipsed (SE) structure has been selected to represent the sandwich family. The 207 PD family is represented by the isomer shown in Figure 2 as it was previously reported 208 to be the most stable of this family at the DFTB level. Two additional structures were 209 considered to account for the T-shaped family, namely the T and Csoa isomers (corre-210 sponding to T4 and Csoa in reference [26]), which only differ by a slight displacement of 211 the top benzene from a symetric position toward a position over a carbon atom. These 212 two structures were previously reported to be degenerated as their energies differ by less 213 than 10 -3 eV at the DFTB level and the present DFT-D3 calculations show a difference 214 of 8 10 -3 eV in favor of the Csoa isomer.215 216

218Table 4 .

 4 the case of the T-shaped family, only one of the two studied isomers could be located on 219 the DFTB potential energy surfaces, namely Csoa with the DFTB-D1 method and T with 220 the DFTB-D2 one. DFT-D3 and ab initio reference calculations agree on the fact that T-221 shaped and PD structures are close in energy and by far more stable than the SE structure. 222 This ordering is not reproduced at the DFTB-D1 nor at the DFTB-D2 level, i.e. for the 223 two dispersion corrections investigated in the absence of atomic charge corrections, 224 because the SE structure is found to be almost degenerated with the PD structure and 225 the T-shaped isomer is found to be the less stable one in both cases. Introducing the 226 WMull charge correction detailed in section 1 with the D1 dispersion (DFTB-D1-WMull) 227 makes the PD structure the most stable but the T-shaped structure remains the less stable. 228 Finally, the DFTB-D2-WMull method provides a correct picture with the T-shaped and 229 PD isomers being close in energy and more stable than the SE isomer. In addition, the 230 binding energies are of the same order as those of the reference calculations. 231 Binding energies for benzene dimer bz 2 in eV. Csoa -0.126 -0.194 -0.192 DFTB-D1-WMull Csoa -0.140 -0.162 -0.148 DFTB-D2 T -0.099 -0.135 -0.132 DFTB-D2-WMull T -0.113 -0.104 -0.086 DFT-D3

3. 4 .Figure 2 .

 42 Figure 2. Benzene dimers bz 2 in vacuum and deposited on graphene

282

  In the present paper, we have reported a new implementation of periodic boundary 283 conditions in the DFTB code deMonNano, as only the Γ-point approximation was 284 available in the previous version of the code. An originality of our scheme is the 285 inclusion of atomic charge corrections which improves the description of intermolecular 286 coulomb interactions. It allows to recover a reasonable description of molecular clusters, 287 as shown in the particular case of benzene dimers in this work. Dispersion corrections are 288 also mandatory for a proper description of such interactions and we have benchmarked 289 two empirical correction schemes. One of them gives the best C-C bond distance in 290 graphene and graphite whereas the second one provides the best interlayer distance in 291 graphite, according to previous reference calculations and experiments, as well as with 292 new DFT calculations performed with the DFT-D3 dispersion corrected functional. 293 Benzene monomer and dimers have been optimized at the DFTB and DFT levels, 294 providing the following similar trends. For the deposition of a single benzene monomer on a graphene sheet, the adsorption of the benzene centered on top of graphene carbon 296 atom or C-C bond leads to almost degenerated structures, by far more stable than 297 the superimposition of the benzene on top of a graphene aromatic cycle. The most 298 stable one at the DFTB level has been selected to build initial conditions for benzene 299 dimers deposition on graphene. The structural energy gap between the most stable 300 T-shaped dimer and less stable Sandwich-like dimer is divided by two when the cluster 301 is supported on graphene. The supported Parallel-Displaced structure appeared to be 302 not stable at the DFTB level, leading to a structure where the two benzene are deposited 303 close to each other on the graphene surface. This structure is the most stable one of 304 our calculations at DFT and DFTB levels, also more stable than the deposition of two 305 benzene monomer at infinite distance, which is not the case of the deposited sandwich or 306 T-shaped dimers. As a conclusion, we have shown the ability of the new implementation 307 to characterize properties of molecular clusters deposited on surfaces.308 Acknowledgments: This work was granted access to the HPC resources of CALMIP supercom-309 puting center under the allocations p18009.

  and the reduction of integrals to one-or two-center

	79
	80	terms :
	81	•	H 0 µ,ν∈α (ρ 0 ) ≈ H 0 µν (ρ α 0 ) ≈ δ µν µα ; the atomic orbital energies of the isolated atom α
	82	•	H 0 µ∈α,ν∈β (ρ 0 ) ≈ H 0 µν (ρ α 0 + ρ	β 0 ) which only depends on the distance between the two
	83		corresponding atomic centers : H 0 µ∈α,ν∈β (r α -r β )

  -5 Hartree are taken into account, in order to limit the When calculations are performed with the WMull scheme, a value of t CH = 0.245 has 115 been determined to provide the atomic charges for the benzene molecule in agreement 116 with reference calculations (see tables and discussion in reference[START_REF] Rapacioli | Correction for dispersion and Coulombic interactions in molecular clusters with density 388 functional derived methods: Application to polycyclic aromatic hydrocarbon clusters[END_REF]).

	113
		number N of boxes involved in the sum.
	117
	118	Regarding convergency criterions, we have used a tolerance of 10 -8 for the atomic
	119	charges during the SCC process and 5.10 -6 Hartree/Bohr for the largest gradient for
		local optimizations.

114 120 2.2. DFT calculations 121 Dispersion corrected DFT calculations were performed under periodic boundary 122 conditions using the Vienna ab initio simulation package (VASP [34-36]) together with 123

  -3 Å . Values of 1.430±0.001 Å and 1.46±0.001 Å were obtained with the DFTB-D1 and DFTB-D2 methods, respectively (see Table1). These 156 values are slighly larger than the value of 1.421 Å previously reported by Zhechkov et al. 157 using the Γ-point approximation and a smaller unit cell [33]. The DFTB-D2 values gives 158 the best agreement with the C-C bond length obtained at the DFT-D3 level (1.425±0.001 159 Å) as well as with the experimental values (1.42 Å).

	160
	161	The graphite bulk has been modeled by including two layers of the previously
	162	defined graphene sheet in the periodic box. In order to determine the appropriate
	163	number of k-points in the z direction (perpendicular to the graphene planes), we have
	164	performed single point energy calculations for an interlayer distance of 3.5 Å, chosen
		because it corresponds to the DFT-D3 one (3.488 Å, see Table

Table 192

 192

	195
	about 35% (∼0.2 eV).

2). The absolute binding energies provided by the DFTB-D2 scheme are in very good 193 agreement with DFT results (appart from LDA) and experimental measurements. The 194 DFTB-D1 scheme gives poorer results, with an overestimation of the binding energies of 196 Regarding the z-separation between the benzene monomer and the graphene sheet (see 197 table 3), DFTB-D1 and DFTB-D2 clearly underestimate the benzene-graphene distance by 198 ∼0.35 Å. However, it should be noted that z-separations calculated with the dispersion 199 corrected DFT functionals also significantly differ from each other by up to 0.25Å.

200

Table 3 .

 3 Z-separation of Benzene on Graphene (in Å).

	Method	a1	a2	a3	a4
	DFTB-D1	3.152 3.133 3.130 3.138
	DFTB-D2	3.080 3.073 3.081 3.054
	DFT-D3	3.465 3.382 3.358 3.415
	ωB97X-D [49]	3.36	3.30	3.35
	vdW-DF1

Table 5 .

 5 Binding energy of benzene dimers bz 2 on graphene in eV. * No value is reported at the DFTB-D2-WMull level for PD as the optimisation led to the a 4cd structure.

	Dissociation ref.		DFTB-D2-WMull			DFT-D3
		T	PD * a 4 -cd	SE	T	PD	a 4 -cd	SE
	graphene + 2bz	-0.586	-0.929 -0.573 -0.605 -0.624 -0.938 -0.583
	graphene + bz 2	-0.473	-	-0.487 -0.467 -0.473	-	-0.478
	graphene@bz + bz -0.135	-0.478 -0.122 -0.156 -0.175 -0.489 -0.134
	graphene@2bz	0.316	-0.027	0.329	0.294	0.275	-0.039	0.316
	The binding energies associated to the optimized structures are reported in

Table 252 5

 252 making use of various choices for the potential energy zeroth reference. In the first 253 line (graphene + 2bz), the reference energy is the sum of the energies for an optimized 254 graphene monolayer and two isolated benzene molecules. It appears that, at the DFT 255 and DFTB levels, the most stable configuration relies on the dissociation of the benzene 256 dimer to form the a 4cd structure. The energetic difference between the T-shaped 257 structure and the less stable SE dimer is twice smaller (0.013 eV vs 0.027 eV at the DFTB 258 level and 0.022 eV vs 0.052 eV at the DFT level) when the dimer is deposited with respect 259 to the gas phase condition. This is probably due to the interaction between the graphene 260 surface and the benzene molecule that is furthest from the surface, which is favored in 261 the sandwich configuration. In the second line (graphene + bz 2 ), the reference energy

262

is that of an isolated graphene sheet plus that of the optimized dimer in its T, PD or SE 263 form, respectively. The gained energy for the non-dissociating dimers (T and SE) are 264 similar (∼0.46-0.49 eV, for DFT and DFTB values), only very slighlty above the binding 265 energy of a single benzene with graphene (0.45 eV for DFT and DFTB values). In the 266 third line (graphene@bz + bz), the reference energy is that of a benzene deposited on 267 a graphene sheet plus that of an isolated benzene. It differs from the isolated dimers 268 by 0.022 eV for the T-shaped structure and 0.036 eV for the SE structure at the DFTB 269 level and 0.018 eV and 0.028 eV at the DFT level. The larger value obtained for the 270 SE dimer can be, again, related to the expected larger interaction energy between the 271 graphene sheet and the second further benzene unit in the SE configuration. The last 272 line (graphene@2bz) compares the binding energies with the one of a system where two 273 benzene molecules would be deposited in a 4 configurations without interaction between 274 them. This configuration appears to be more stable than those corresponding to the 275 deposition of a T-shaped or SE dimer. The negative sign obtained for the a 4cd structure 276

Version February 10, 2022 submitted to Journal Not Specified https://www.mdpi.com/journal/notspecified

Conflicts of Interest:

The authors declare no conflict of interest.

Physical


Laurence Jacquet
email: laurence.jacquet@u-cergy.fr.

Etienne Lehmann
email: etienne.lehmann@gmail.com.

Pierre Boyer

Vidar Christiansen

Guy Laroque

Emmanuel Saez

Laurent Simula

Stefanie Stantcheva

Kevin Spiritus

Alain Trannoy

Optimal tax problems with multidimensional heterogeneity: A mechanism design approach *

Keywords: Optimal taxation, mechanism design, multidimensional screening problems, allocation perturbation

We propose a new method, that we call an allocation perturbation, to derive the optimal nonlinear income tax schedules with multidimensional individual characteristics on which taxes cannot be conditioned. It is well established that, when individuals differ in terms of preferences on top of their skills, optimal marginal tax rates can be negative. In contrast, we show that with heterogeneous behavioral responses and skills, one has optimal positive marginal tax rates, under utilitarian preferences and maximin.

I Introduction

Since [START_REF] Mirrlees | An Exploration in the Theory of Optimum Income Taxation[END_REF], the mechanism design method has been largely used to obtain optimal tax formulas when individuals differ along a single dimension of heterogeneity, typically skills, and when income taxation cannot be conditioned on this dimension of heterogeneity. 1In this paper, we extend [START_REF] Mirrlees | An Exploration in the Theory of Optimum Income Taxation[END_REF]'s mechanism design method to a situation where taxpayers differ along many characteristics that are unobserved or on which taxation cannot be conditioned. 2 For instance, taxpayers with the same income level can differ in terms of unobserved preferences for work, gender, number and age of the children, cultural background, geographical area where one lives, ethnicity or age. Taxpayers can then exhibit heterogeneous labor supply elasticities in addition to distinct skill levels.

In this paper, we develop a new method, that we call "allocation perturbation", to solve multidimensional screening problems. We show that individuals of different groups pooled at the same income level are characterized by the same marginal rate of substitution between pre-tax and after-tax income. Intuitively, individuals of distinct groups who earn the same income level face the same marginal tax rate. From the individual maximization program, we know that identical marginal tax rates imply identical marginal rates of substitution. Using this equality in marginal rates of substitution together with a single-crossing condition within each group, we can fully characterize an incentive-compatible allocation. We then apply variational calculus to small perturbations of this allocation and derive the necessary conditions of the optimal allocation. To do so, we take the Gâteaux derivatives of the associated Lagrangian and equate them to zero. Rearranging terms of these conditions gives the optimal tax formula.

In several recent papers, tax perturbation methods are proposed to solve multidimensional screening problems (e.g., [START_REF] Hendren | Measuring Economic Efficiency Using Inverse-Optimum Weights[END_REF] and [START_REF] Sachs | Nonlinear Tax Incidence and Optimal Taxation in General Equilibrium[END_REF] and [START_REF] Jacquet | Optimal redistributive taxation with both extensive and intensive responses[END_REF] for the most recent). In contrast to our approach, all these papers are built upon the tax perturbation approach initiated by [START_REF] Piketty | La Redistribution fiscale contre le ch ômage[END_REF] and [START_REF] Saez | Using Elasticities to Derive Optimal Income Tax Rates[END_REF] which consists in computing all responses to small tax reforms, to sum them up and equate them to zero in order to obtain the optimal tax schedule.

Tax perturbation approaches easily allow one to obtain optimal tax formulas in terms of meaningful sufficient statistics. In contrast, we derive the optimal income tax schedule in terms of policy invariant functions (i.e. skill density and the derivatives of the individual and social utility functions) and rewrite it in terms of sufficient statistics. Just like the tax formula in terms of sufficient statistics, our formula in terms of policy invariant functions is not closedform because it depends on the allocation where these functions are evaluated. We use it to sign the optimal marginal tax rates. With multidimensional heterogeneity, the literature has highlighted that negative marginal tax rates can be optimal. In [START_REF] Cuff | Optimality of workfare with heterogeneous preferences[END_REF], Boadway et al. (2002b), [START_REF] Brett | How optimal nonlinear income taxes change when the distribution of the population changes[END_REF], [START_REF] Choné | Negative Marginal Tax Rates and Heterogeneity[END_REF] and [START_REF] Lockwood | De Gustibus non est Taxandum: Heterogeneity in Preferences and Optimal Redistribution[END_REF], individuals differ along their skills and preferences for leisure, and the social planner has weighted utilitarian preferences. In this context, individuals who pool at the same income level are characterized by different marginal social welfare weights. Therefore, the mean social welfare weight may not be decreasing with income as it is the case with one dimension of heterogeneity and e.g. utilitarian preferences. This opens the way to negative marginal tax rates. One may wonder whether heterogeneous labor supply elasticities for individuals who earn the same level of income could also open the way to negative marginal tax rates. We show that, under utilitarian preferences, marginal tax rates are positive. This result also prevails under maximin.

Our method can easily be extended to include participation decisions (Saez, 2002[START_REF] Kleven | The Optimal Income Taxation of Couples[END_REF][START_REF] Jacquet | Optimal redistributive taxation with both extensive and intensive responses[END_REF], migration decisions [START_REF] Lehmann | Tax Me if You Can! Otimal Nonlinear Income Tax Between Competing Governments[END_REF][START_REF] Blumkin | International tax competition: zero tax rate at the top re-established[END_REF] or sectoral decisions [START_REF] Rothschild | Redistributive Taxation in the Roy Model[END_REF][START_REF] Scheuer | Adverse selection in credit markets and regressive profit taxation[END_REF][START_REF] Gomes | Differential Taxation and Occupational Choice[END_REF], simply following [START_REF] Jacquet | Optimal redistributive taxation with both extensive and intensive responses[END_REF]. The tax formulas then simply incorporate new terms implied by the participation margin.

The paper is organized as follows. In Section II, we present the model. In Section III, we characterize any incentive compatible allocation thanks to a monotonicity constraint, a differential equation and a pooling function. In Section IV, we detail the allocation perturbation method and derive the optimal tax schedule. In Section V, we give a sufficient condition for optimal marginal tax rates to be positive. Section VI concludes.

II Model

II.1 Individuals

Taxpayers differ along their skill level w ∈ R * + and along some characteristics denoted θ ∈ Θ. We call a group a subset of individuals with the same θ. 3 We assume that the set of groups Θ is measurable with a cumulative distribution function (CDF) denoted µ(•). The set Θ can be finite or infinite and may be of any dimension and is compact. The distribution µ(.) of the population across the different groups may be continuous, but it may also exhibit mass points. Among individuals of the same group θ, skills are distributed according to the conditional skill density f (•|θ) which is positive and differentiable over the support R * + . The conditional CDF is denoted F(w|θ) def ≡ w 0 f (x|θ)dx. We do not make any restriction on the correlation between w or θ. We normalize to unity the total size of the population.

Every worker derives utility from consumption c ∈ R + and disutility from effort. Effort captures the quantity as well as the intensity of labor supply. More effort implies higher pretax income y ∈ R + (for short, income hereafter). The government levies a tax T(.) which depends on income y only. Consumption c is related to income y through the tax function T(y) according to c = y -T(y).

Assumption 1. The utility function is additively separable and takes the form:

U (c, y; w, θ) = u(c) -υ(y; w, θ) with : u , υ y , υ y,y > 0 > υ w , u ≤ 0
The convexity of the indifference curve is ensured by assuming that υ y,y > 0 ≥ u . Additive separability is standardly used in the adverse selection literature with multidimensional heterogeneity, e.g., [START_REF] Rochet | The taxation principle and multi-time Hamilton-Jacobi equations[END_REF], [START_REF] Wilson | Nonlinear Pricing[END_REF], Rochet and Choné (1998), [START_REF] Rochet | Nonlinear Pricing with Random Participation[END_REF]. As in the model with one-dimensional heterogeneity [START_REF] Mirrlees | An Exploration in the Theory of Optimum Income Taxation[END_REF], this assumption is necessary to sign optimal marginal tax rates (Proposition 3). Note that additive separability can be seen as a rather restrictive assumption, in a multidimensional context, since the group θ matters only for the marginal disutility of income but not for the marginal utility of consumption.

A worker of type (w, θ), facing y → T(y), solves:

U(w, θ) def ≡ max y u(y -T(y)) -υ(y; w, θ) (1)
We call Y(w, θ) the solution to program (1), C(w, θ) = Y(w, θ) -T(Y(w, θ)) the consumption of a worker of type (w, θ) and U(w, θ) her gross utility (or maximized level of utility). When the tax function is differentiable, the first-order condition associated to (1) implies that:

1 -T (Y (w, θ)) = M (C (w, θ) , Y (w, θ) ; w, θ) (2)
where:

M (c, y; w, θ) def ≡ υ y (y; w, θ) u (c) (3)
denotes the marginal rate of substitution between income and consumption. For a worker of a given type, the left-hand side of Equation (2) corresponds to the marginal gain of income after taxation while the right-hand side corresponds to the marginal cost of income in monetary terms.

We impose the single-crossing (Spence-Mirrlees) condition that, within each group of workers endowed with the same θ, the marginal rate of substitution is a decreasing function of the skill level, i.e. that higher-skilled workers find it less costly to increase their income y:

Assumption 2 (Within-group single-crossing condition). For each θ ∈ Θ, and each (c, y) ∈

R + × R + , function w → M (c, y; w, θ) is differentiable with ∀w ∈ R * + , M w < 0 i.e. υ y,w < 0, lim w →0 M (c, y; w, θ) = +∞ and lim w →∞ M (c, y; w, θ) = 0, i.e. lim w →0 υ y,w (y; w, θ) = +∞ and lim w →∞ υ y,w (y; w, θ) = 0.

II.2 Government

The government's problem consists in finding the tax schedule T(•) that maximizes a social welfare function,

θ∈Θ,w∈R * + Φ (U(w, θ); w, θ) f (w|θ)dw dµ(θ), (4)
where Φ(U; w, θ) is an increasing transformation of individuals' utility levels U. This welfarist specification allows Φ to vary with type (w, θ) which makes it very general. Weighted utilitarian preferences (as in Weymark (1987)) are obtained with Φ(U; w, θ) ≡ ϕ(w, θ) • U with weights ϕ(w, θ) depending on individual characteristics. The objective is utilitarist if ϕ(w, θ) is constant and Φ(U; w, θ) ≡ U and it turns out to be maximin (or Rawlsian) if ϕ(w, θ) = 0 ∀w > 0.

When Φ(U; w, θ) does not vary with its two last arguments and is concave in individual utility (Φ UU ≤ 0), we obtain a Bergson-Samuelson criterion which is a concave transformation of utility.

When maximizing social welfare, the government takes into account the individual labor supply decisions (1) and the budget constraint,

θ∈Θ,w∈R * + T (Y(w, θ)) f (w|θ)dw dµ(θ) = E, (5)
where E ≥ 0 is an exogenous amount of public expenditures. Let λ > 0 denote the shadow price of public funds.

III Incentive constraints and pooling function

According to the taxation principle [START_REF] Hammond | Straightforward Individual Incentive Compatibility in Large Economies[END_REF][START_REF] Guesnerie | A Contribution to the Pure Theory of Taxation[END_REF], it is equivalent for the government to select a nonlinear tax schedule taking into account labor supply decisions described in (1), or to directly select an allocation (w, θ) → (C(w, θ), Y(w, θ)) that verifies the incentive constraints:

∀w, θ, w , θ ∈ (R * + × Θ) 2 U (C(w, θ), Y(w, θ); w, θ) ≥ U C(w , θ ), Y(w , θ ); w, θ . (6)
The incentive constraints (6) impose that workers of type (w, θ) prefer the bundle designed for them rather than the bundles (C(w , θ ), Y(w , θ )) designed for individuals of any other type (w , θ ).

With multidimensional heterogeneity, incentive constraints (6) can be rewritten as incentive constraints within each group θ and incentive constraints across distinct groups. In a first subsection, we show that the within-group incentive constraints can be reformulated as a monotonicity constraint and a differential equation, as in the model with one dimension of heterogeneity. In a second subsection, we show that, across groups, to guarantee incentive compatibility, taxpayers from distinct θ-groups who earn the same level of income need to face the same marginal rate of substitution between consumption and income. To satisfy this equalization of marginal rates of substitution across groups, we define a pooling function that describes the level of skill required in each θ group to obtain a given level of income. We show that we have an incentive-compatible allocation under the additional condition that, for each θ, the pooling function be smoothly increasing in skill.

III.1 Incentive constraints within groups

An incentive-compatible allocation has to satisfy (6). It thus has to verify for each group θ the following set of "within-group incentive constraints".

∀(w, w, θ) ∈ R 2 + × Θ U (C(w, θ), Y(w, θ); w, θ) ≥ U (C( w, θ), Y( w, θ); w, θ) . ( 7)

Under the within-group single-crossing assumption 2, the set of within-group incentive constraints can be transformed into a monotonicity constraint and a differential equation that we

give in Lemmas 1 and 2 below.

Lemma 1. Under Assumption 2, the function w → Y(w, θ) is nondecreasing for each θ ∈ Θ.

Assumption 2 (i.e. the within-group single-crossing assumption) implies that, in the same group θ, the indifference curves of workers with a lower skill are steeper than the ones of workers with a higher level of skill, as in one dimensional tax models. We therefore skip the proof of Lemma 1 which is a simple reformulation of the usual proof in the one dimensional framework. Note that Y(•; θ) being nondecreasing, it may exhibit discontinuities over a countable set and it may also exhibit bunching if individuals in the same group but endowed with different skill levels earn the same income. We follow a first-order approach and consider only smooth allocations where these two pathologies do not arise.To do so, we make a smoothness assumption. As a preamble, we define smoothly increasing functions 4 and then give the smoothness assumption.

Definition 1. A function a : R + → R is "smoothly increasing"if it is differentiable with ∀x ∈ R + , a (x) > 0, a (0) = 0 and lim

x →∞ a (x) = +∞.

Assumption 3 (Smooth allocations). In each group θ, w → Y(w, θ) is smoothly increasing.

According to Assumption 3,5 for each income level y ∈ R + and for each group θ ∈ Θ, there exists a single skill level w such that only individuals of that skill level within group θ earn income y = Y(w, θ). The following lemma provides the first-order incentive constraints within group θ reformulated as a differential equation, where the dot above a variable stands for the partial derivative of this variable with respect to skill w.

Lemma 2. Under Assumptions 1 and 3, for each θ, the mapping w → U(w, θ) is differentiable with:

U(w, θ) = -v w (Y(w, θ); w, θ) . (8a)
Moreover, Equation (8a) is equivalent to:

Ċ(w, θ) Ẏ(w, θ) = M (C(w, θ), Y(w, θ); w, θ) . ( 8b
)
We skip the proof since it simply consists in adapting that of the one-dimensional tax models, which consists in applying the envelope theorem to Equation (1), see e.g., [START_REF] Salanié | The Economics of Taxation[END_REF].

Integrating (8a) leads to:

U(w, θ) = U(0, θ) - w 0 v w (Y(x, θ); x, θ) dx. (8c)
The "first-order" approach we follow is usual in one-dimensional tax models (see e.g. Salanié ( 2011)), since it takes into account the first-order incentive-compatibility constraints (8a) and assumes the second-order incentive-compatibility constraints are satisfied (Assumption 3) as standardly observed on data.6 .

III.2 Incentive constraints across groups and pooling function

We now describe how the various within-group allocations ω → (Y(ω, θ), C(ω, θ)) need to be set to be mutually incentive-compatible and to verify the full set of incentive constraints (6). This is the pooling issue that we now address.

Choose a reference group θ 0 ∈ Θ, a skill level w and another group θ. Individuals of type (w, θ 0 ) earn income Y(w, θ 0 ). According to the smoothness assumption 3, each group-specific

allocation Y(•, θ) : w → Y(w, θ
) is an increasing one-to-one function that maps the positive real line onto itself. Therefore, there must exist a single skill level, hereafter denoted W(w, θ), so that individuals of the other group θ endowed with that skill level W(w, θ) must get the same income level Y(w, θ 0 ) as individuals of type (w, θ 0 ), i.e. Y(W(w, θ), θ) = Y(w, θ 0 ). We call W(., .) the pooling function. For each θ ∈ Θ, the pooling function combines two smoothly increasing functions, namely ω

Y(•,θ 0 ) -→ Y(ω, θ 0 ) Y -1 (•,θ) -→ W(ω, θ).
The pooling function is therefore also a smoothly increasing function in skill w. It obviously verifies W(w, θ 0 ) ≡ w. Provided that the allocation is incentive-compatible, it is not possible from (6) that individuals of type (W(w, θ), θ) and individuals of type (w, θ 0 ) obtain the same income Y(w, θ 0 ) but distinct consumption levels. Therefore, for each (w, θ), we must simultaneously have

Y(W(w, θ), θ) ≡ Y(w, θ 0 ) and C(W(w, θ), θ) ≡ C(w, θ 0 ). (9)
Individuals who earn the same income level face the same marginal tax rate, according to Equation ( 2), hence have the same consumption level. The skill levels of individuals who earn a given income level are implicitly determined by the equality of their marginal rates of substitution which is our pooling condition.

Lemma 3. Under Assumptions 2 and 3, along an incentive-compatible allocation, the bundle designed for individuals of type (W(w, θ), θ) coincides with the bundle (C(w, θ 0 ), Y(w, θ 0 )) designed for individuals of type (w, θ 0 ), where W(w, θ) verifies the following pooling condition:

M (C(w, θ 0 ), Y(w, θ 0 ); w, θ 0 ) = M (C(w, θ 0 ), Y(w, θ 0 ); W(w, θ), θ) . ( 10
) Proof According to Assumption 2, M (C(w, θ 0 ), Y(w, θ 0 ); w, θ 0 ) = M (C(w, θ 0 ), Y(w, θ 0 ); ω, θ)
admits exactly one solution in ω. Differentiating in w both sides of equalities in (9) and rearranging terms leads to:

Ċ(W(w, θ), θ) Ẏ(W(w, θ), θ) = Ċ(w, θ 0 ) Ẏ(w, θ 0 ) .
According to Lemma 2, Equation (8b) holds, which implies (10).

One can retrieve the entire incentive-compatible allocation for all groups if one knows the pooling function W(•, •) and the allocation ω → (Y(ω, θ 0 ), C(ω, θ 0 )) designed for the reference group.

Thanks to the pooling condition, in the following lemma, we provide a sufficient condition for the allocation to be incentive-compatible. The proof is in Appendix A.1.

Lemma 4. Under Assumption (2), let w → (C(w, θ 0 ), Y(w, θ 0 )) be a within-group allocation that verifies Assumption 3 and the within-group incentive-compatible Equation (8b). For each w ∈ R + and each group θ ∈ Θ, let W(w, θ) be the unique skill level ω that solves the pooling condition

M (C(w, θ 0 ), Y(w, θ 0 ); w, θ 0 ) = M (C(w, θ 0 ), Y(w, θ 0 ); ω, θ).
There exists a unique incentive-compatible allocation (w, θ) → (C(w, θ), Y(w, θ)) the restriction of which to group θ 0 is w → (C(w, θ 0 ), Y(w, θ 0 ))

and it verifies Assumption 3 if and only if, for each θ, w → W(w, θ) is smoothly increasing.

Lemma 4 guarantees that if w → Y(w, θ) is smoothly increasing in w and if, for each θ, the pooling function denoted W(w, θ) is also smoothly increasing in w, then the allocation is incentivecompatible. Assumption 3 together with the assumption that W(•, θ) is smoothly increasing plays here a role similar to the assumption that the second-order incentive compatibility condition is satisfied with one dimension of heterogeneity. In what follows, we therefore select the allocation only for the reference group θ 0 and assume that the triggered allocations for the other groups verify Assumption 3. From Assumption 1 and Equation (3), the pooling condition (10) can be rewritten as:

v y (Y(w, θ 0 ); w, θ 0 ) = v y (Y(w, θ 0 ); W(w, θ), θ) . ( 11
)
The pooling function W(•, θ) that enables to retrieve (C(•, θ), Y(•, θ)) from the allocation of the reference group (C(•, θ 0 ), Y(•, θ 0 )) depends on Y(•, θ) and not on C(•, θ) thanks to the additive separability of the utility function. This is necessary to apply our allocation perturbation method. Interestingly, the pooling function is endogenous so that individuals with the same income level can have distinct behavioral elasticities. This is a major difference with the previous literature where characteristics are aggregated along a single dimension (which can encapsulate general equilibrium effects on the wage distribution for instance), e.g. Boadway et al. (2002a), [START_REF] Brett | How optimal nonlinear income taxes change when the distribution of the population changes[END_REF], [START_REF] Rothschild | Redistributive Taxation in the Roy Model[END_REF]Scheuer (2013, 2016), Scheuer (2013) and [START_REF] Gomes | Differential Taxation and Occupational Choice[END_REF]. Such an aggregation indeed prevents individuals who earn the same income to differ in terms of behavioral elasticities.

IV Allocation perturbation and optimal tax schedule

In this section, in order to obtain the necessary conditions for the optimal tax schedule, we propose a new method that we call an "allocation perturbation". We first reformulate the government problem taking into account the pooling function we introduced in the previous section. We second motivate and summarize our methodology. We then present the optimal tax schedule and the technical details behind its derivation. Eventually, we rewrite our optimal tax formula in terms of sufficient statistics.

Let C ( û, y; w, θ) denote the consumption level the government needs to provide to a worker of type (w, θ) who earns y to ensure she enjoys a given utility level û. The function C (•, y; w, θ)

is the reciprocal of U (•, y; w, θ) and:

C u ( û, y; w, θ) = 1 u (c) and C y ( û, y; w, θ) = v y (y; w, θ) u (c) (12)
where the various derivatives are evaluated at c = C ( û, y; w, θ). Define the Lagrangian associated to the planner's problem as

L def ≡ Y(w, θ) -C (U(w, θ), Y(w, θ); w, θ) + Φ (U(w, θ); w, θ) λ f (w|θ)dw dµ(θ). ( 13
)
where all terms are expressed in monetary units.

Under Assumption 3, thanks to Lemma 4, the government then simply chooses, among the set of smooth allocations, the best one that verifies the first-order incentive constraint (8a) for each group, and the pooling condition (10). When the pooling function associated with the solution is, for each group θ, smoothly increasing in skill w, the found solution also solves the problem with all incentive constraints.

To solve this type of problem with one-dimensional unobserved heterogeneity, one typically constructs a Hamiltonian and one applies the Pontryagin principle. In our multidimensional environment, the pooling condition (10) induces constraints on state and control variables which hold at endogenous skill levels. In this context, we rather propose using the calculus of variation and consider a set of perturbations of the allocation in the reference group.

The cornerstone of our method is the pooling condition (10) that we use to deduce how the allocations in the other groups are perturbed.

IV.1 Allocation perturbation: Methodology

Consider a perturbation of the allocation for a given "reference" group that preserves incentive compatibility. In this group, change all incomes that correspond to a small interval of skills (wδ, w) and do not modify incomes outside this interval. More precisely, perturb incomes in

the reference group Y(x, θ 0 ) to obtain Y(x, θ 0 ) + t ∆ Y (x, θ 0 ; δ) where ∆ Y (•, θ 0 ; δ) is a continuously differentiable function defined on R + such that ∆ Y (•, θ 0 ; δ) is positive for x ∈ (w -δ, w) and
is nil otherwise, and where t ∈ R is an algebraic magnitude. Thanks to the pooling function and Lemma 4, one can deduce the impact of this perturbation in all groups (within and outside the skills interval where the perturbation takes place). Take the Gâteaux derivative of the perturbed Lagrangian associated to the government's optimization problem with respect to t, at t = 0, and equate it to zero. One obtains an equation that characterizes the optimal marginal tax rates. The other necessary condition of the optimal tax schedule is obtained by perturbating the government's Lagrangian as follows. Uniformly modify utility (without modifying incomes) of all taxpayers by an amount ∆. Then take the Gâteaux derivative of the perturbed Lagrangian with respect to ∆, at ∆ = 0, and equate it to zero.

IV.2 Derivation of the optimal tax schedule

To save on notations, we from now on use the more compact notation [x, θ] when the vari- ous functions are evaluated for types (x, θ) at income Y(x, θ), utility U(x, θ) and consumption C(x, θ). The optimality condition with multidimensional heterogeneity is presented in the following proposition.

Proposition 1. Under Assumptions 1, 2 and 3, the optimal tax schedule verifies:

T (Y(w, θ 0 )) 1 -T (Y(w, θ 0 )) θ∈Θ v y [W(w, θ), θ] -W(w, θ) v yw [W(w, θ), θ] W(w, θ) f (W(y, θ)|θ) dµ(θ) (14a) = u (C(w, θ 0 )) θ∈Θ,x≥W(w,θ) 1 u (C(x, θ)) - Φ U (U(x, θ); x, θ) λ f (x|θ)dx dµ(θ)
for all skill w and θ∈Θ,w∈R

+ Φ U (U(w, θ); w, θ) λ - 1 u (C(w, θ)) f (w|θ)dw dµ(θ) = 0. ( 14b
)
Proof (Allocation perturbation) To derive (14a) at a given skill level w, we consider a set of allocation perturbations, indexed by t ∈ R and δ ∈ R + , that we denote Ĉ(w, θ; t, δ) , Ŷ(w, θ; t, δ) and Û(w, θ; t, δ) def ≡ U Ĉ(w, θ; t, δ), Ŷ(w, θ; t, δ); w, θ ) where t stands for the size of the perturbation, and δ is the length of the skill interval where, in the reference group, the perturbation of incomes takes place. Following Lemma 4, we define the allocation perturbations from their restriction to the reference group θ 0 and then study the impact of these perturbations on the allocation in every other group. The perturbations of incomes in the reference group are defined by:

Ŷ(x, θ 0 ; t, δ) def ≡ Y(x, θ 0 ) + t ∆ Y (x, θ 0 ; δ)
where ∆ Y (•, θ 0 ; δ) is a continuously differentiable function defined on R + such that ∆ Y (•, θ 0 ; δ) > 0 for x ∈ (wδ, w) and ∆ Y (•, θ 0 ; δ) = 0 otherwise. Incomes in the reference group remain unchanged outside the skill interval (wδ, w) and are increased (decreased) inside the skill interval (wδ, w) if t > 0 (if t < 0), as illustrated in Figure 1. It is worth noting that the perturbed income function remains differentiable with respect to skill w since ∆ Y (•, •, δ) is differentiable. Moreover, from Assumption 3, Y(•, θ 0 ) admits a positive derivative everywhere, so Ẏ(•, θ 0 ) is bounded away from 0 for all x ∈ [wδ, w]. Therefore, provided that t is small enough, which

x Y(x, θ 0 ) w -δ w Ŷ(x, θ 0 ; t, δ) def ≡ Y(x, θ 0 ) +t∆ Y (x, θ 0 ; δ)

Initial allocation Perturbed allocation

Figure 1: The perturbation of incomes in the reference group θ 0 we assume in the rest of the proof, Ŷ(•, θ 0 ; t, δ) has also a positive derivative everywhere and therefore verifies Assumption 3.

Let us in addition assume that the utility of the lowest skilled individuals in the reference group U(0, θ 0 ; t, δ) is not perturbed and write it as U(0, θ 0 ). Therefore, according to (8c), the perturbed utility function in the reference group is

Û(x, θ 0 ; t, δ) def ≡ U(0, θ 0 ) - x 0 υ w Ŷ(ω, θ 0 ; t, δ); ω, θ 0 dω. ( 15a
)
From the pooling condition (10), as incomes Y(•, θ 0 ; t, δ) in the reference group remain unchanged outside the skill interval (wδ, w), the pooling function W(•, θ 0 ; t, δ) is not perturbed outside the skill interval (wδ, w). Therefore, incomes Y(•, θ; t, δ) in any group θ are not modified outside the skill interval (W(wδ, θ), W(w, θ)), and we must have (See Figure 2)

Ŷ(x, θ; t, δ) = Y(x, θ) if x ∈ [0, W(w -δ, θ)] ∪ [W(w, θ), +∞) . (15b) x Y(x, θ) W(w -δ, θ) W(w, θ)
For t > 0

Initial allocation Perturbed allocation

Figure 2: The perturbation of incomes in the other groups Since incomes in the reference group are not perturbed for all skill x below wδ, the pooling function is also unchanged below wδ, so that the same types remain pooled together. Hence we get in group θ that for all x ≤ W(wδ, θ):

Ĉ(x, θ; t, δ) = C(x, θ) and Û(x, θ; t, δ) = U(x, θ). ( 15c
)
For all skills x > W(wδ, θ), the change in utility obtained using the first-order incentive constraint (8c) is:

Û(x, θ; t, δ) -U(x, θ) = - x 0 υ w Ŷ(ω, θ; t, δ); ω, θ -υ w (Y(ω, θ); ω, θ) dω. (15d)
Since incomes Ŷ(•, θ; t, δ) are only perturbed inside (W(wδ, θ)), W(w, θ)), for all skills x that belong to this interval, using (15b), we get:

Û(x, θ; t, δ) -U(x, θ) = x W(w-δ,θ) υ w (Y(ω, θ); ω, θ) -υ w Ŷ(ω, θ; t, δ); ω, θ dω. (15e)
Moreover, for all skills x above W(w, θ), we have:

Û(x, θ; t, δ) -U(x, θ) = W(w,θ) W(w-δ,θ) υ w (Y(ω, θ); ω, θ) -υ w Ŷ(ω, θ; t, δ); ω, θ dω. (15f)
Hence utility in the other group does not change below W(wδ, θ) and changes by a uni-

w U(x, θ) W(w -δ, θ) W(w, θ)
For t > 0

Initial allocation Perturbed allocation

Figure 3: The perturbation of utilities form amount above W(w, θ), as illustrated in Figure 3. As incomes above skill W(w, θ) are unchanged, this implies that, for all skill x above W(w, θ), the modifications in utility U(x, θ) occur only through changes of the utility u(C(x, θ)) derived from consumption. Using (15f), this utility therefore changes uniformly by:

u Ĉ(x, θ; t, δ) -u (C(x, θ)) = W(w,θ) W(w-δ,θ) υ w (Y(ω, θ); ω, θ) -υ w Ŷ(ω, θ; t, δ); ω, θ dω (15g)
which determines the perturbation of consumption for skill levels above W(w, θ). We now determine how the perturbations of incomes Y(•, θ) in each group within the skill interval (W(wδ, θ), W(w, θ)) need to be set to ensure that the perturbed allocations remain incentivecompatible. For that purpose, we note that for all skill levels x above w, as incomes in the reference group are not perturbed, the pooling function is also unchanged, so that the same types remain pooled together. Hence, according to (9):

∀t, ∀x ≥ w Ŷ(W(x, θ), θ; t, δ) = Ŷ(x, θ 0 ; t, δ) and Ĉ(W(x, θ), θ; t, δ) = Ĉ(x, θ 0 ; t, δ).
This implies that, in all groups, the uniform change in utility that occurs for all skill levels above W(w, θ) must be identical across groups, so that: u Ĉ(x, θ 0 ; t, δ)u (C(x, θ 0 )) = u Ĉ(W(x, θ), θ; t, δ)u (C(W(x, θ), θ)), and so, using (8c) and (15g), we obtain:

w w-δ υ w Ŷ(ω, θ 0 ; t, δ); ω, θ -υ w (Y(ω, θ 0 ); ω, θ) dω (15h) = W(w,θ) W(w-δ,θ)
υ w Ŷ(ω, θ; t, δ); ω, θυ w (Y(ω, θ); ω, θ) dω.

The latter equation links the perturbed incomes Ŷ(•, θ; t, δ) in all groups within the interval of skills (W(wδ, θ), W(w, θ)) and the perturbed incomes Ŷ(•, θ 0 ; t, δ) in the reference group.

The perturbed Lagrangian is:

L (t, δ) def ≡ θ∈Θ,w∈R + Ŷ(w, θ; t, δ) -C Û(w, θ; t, δ), Ŷ(w, θ; t, δ); w, θ (16) + Φ Û(w, θ; t, δ); w, θ λ f (w|θ)dw dµ(θ).
If the allocation is optimal, the derivative of this Lagrangian with respect to t must be nil at t = 0. In Appendix A.2, we show that the limit of ∂ L /∂t when δ goes to zero leads to:

θ∈Θ 1 - υ y (Y(W(w, θ), θ); W(w, θ), θ) u (C(W(w, θ), θ) υ yw (Y(W(w, θ), θ); W(w, θ), θ) f (W(w, θ)|θ) dµ(θ) (17) = θ∈Θ,x≥W(w,θ) Φ U [x, θ] λ - 1 u [x, θ] f (x|θ)dx dµ(θ).
Using ( 2), (3), Y(W(w, θ), θ) ≡ Y(w, θ 0 ) and C(W(w, θ), θ) ≡ C(w, θ 0 ), we can rewrite (17) as:

T (Y(w, θ 0 )) θ∈Θ 1 υ yw (Y(W(w, θ), θ); W(w, θ), θ) f (W(w, θ)|θ) dµ(θ) = θ∈Θ,x≥W(w,θ) Φ U [x, θ] λ - 1 u [x, θ] f (x|θ)dx dµ(θ).
Using again (2), ( 3) and ( 9) leads to (14a).

We now derive (14b). Consider a set of allocation perturbations indexed by ∆ ∈ R and denoted ( C(w, θ; ∆), Ỹ(w, θ; ∆), Ũ(w, θ; ∆) def ≡ U C(w, θ; ∆), Ỹ(w, θ; ∆); w, θ ), which consist, for each type (x, θ) ∈ R + × Θ, in no change in Y(x, θ) and in a uniform change in U(x, θ), therefore in u(C(x, θ)) by an amount ∆. Hence, we get for each ∆ that Ũ(w, θ; ∆) def ≡ U(w, θ) + ∆, Ỹ(w, θ; ∆) def ≡ Y(w, θ) and C(w, θ; ∆) def ≡ C Ũ(w, θ; ∆), Ỹ(w, θ; ∆); w, θ . These perturbations preserve incentive-compatibility (6). According to (13), the perturbed Lagrangian can be written as

L (∆) def ≡ Ỹ(w, θ; ∆) -C Ũ(w, θ; ∆), Ỹ(w, θ; ∆); w, θ + Φ Ũ(w, θ; ∆); w, θ λ f (w|θ)dw dµ(θ).
If the allocation is optimal, the above perturbations do not affect the Lagrangian. Thus, by equating the Gâteaux derivative of the Lagrangian in the direction described by the above perturbations (i.e. the derivative of the perturbed Lagrangian L (•) with respect to ∆, at ∆ = 0) to zero, we obtain an equation that characterizes the optimal tax system. Using the first equality in ( 12), this Gâteaux derivative of the Lagrangian is

L (0) = θ∈Θ,x∈R + Φ U (U(x, θ); x, θ) λ - 1 u (C(x, θ)) f (x|θ)dx dµ(θ).
Equating this derivative to zero leads to (14b).

The tax formula (14a) of Proposition 1 generalizes the [START_REF] Mirrlees | An Exploration in the Theory of Optimum Income Taxation[END_REF] formula and ABC formula described in [START_REF] Diamond | Optimal Income Taxation: An Example with U-Shaped Pattern of Optimal Marginal Tax Rates[END_REF] to multidimensional individual characteristics with an efficiency, an equity and a distribution term. Using the Hamiltonian method, we show in Appendix A.3 that, under a single dimension of heterogeneity w, Equations ( 14a) and (14b) simplify to (18a)-(18b).

Lemma 5. When the unobserved heterogeneity has only one dimension, the optimal tax schedule satisfies:

T (Y(w)) 1 -T (Y(w)) • v y [w] -w v yw [w] w f (w) = u [w] ∞ w 1 u [x] - Φ U [x] λ f (x) dx (18a) 0 = ∞ 0 1 u [x] - Φ U [x] λ f (x) dx. ( 18b
)
Comparing these equations with Equations ( 14a) and ( 14b) makes clear that reducing the tax problem to one dimension of heterogeneity implies that the integrals over θ-groups disappear. With multidimensional heterogeneity, one needs to aggregate the terms of the formula for individuals of the different groups who pool at the same level of income. This is made possible thanks to our characterization of the pooling function in Lemmas 3 and 4.

Comparing Proposition 1 with Lemma 5 makes it clear that, with multidimensional heterogeneity, the composition of the population at each income level plays a role through the weight

W(y, θ) f (W(y, θ)|θ) dµ(θ)
. This weight is not simply the density of people in each group θ who earn the relevant level of income but the product of the skill levels W(y, θ) times the corresponding densities f (W(y, θ)|θ) dµ(θ). We know that taxpayers who earn the same income come from distinct groups θ have distinct skill levels W(y, θ). Moreover, at each income level, the composition of the population (e.g. the proportion of men versus women) usually varies between the actual economy where it is estimated and the optimal economy. Therefore, the weights W(y, θ) f (W(y, θ)|θ) dµ(θ) take distinct values in the actual and optimal economies. In [START_REF] Jacquet | Optimal redistributive taxation with both extensive and intensive responses[END_REF], we show that ignoring these so-called composition effects significantly biases the optimal income tax schedule.

IV.3 Rewriting the optimal tax schedule in terms of sufficient statistics

In this section, we rewrite the optimal tax formula of Proposition 1 in terms of meaningful sufficient statistics. To obtain the latter, we first define a set of individual elasticities. For this purpose, we assume that the tax function T(•) is twice differentiable which ensures that the first-order condition associated to the individual maximization program (1) is differentiable.

We also assume that, for all (w, θ) ∈ R * + × Θ, the second-order condition associated to the individual maximization program holds strictly which guarantees it is invertible in y. This allows one to apply the implicit function theorem to the first-order condition associated to the individual maximization program. We also assume that, for all (w, θ) ∈ R * + × Θ, the function y → U (y -T(y), y; w, θ) admits a unique global maximum over R + . This ensures that after a change in the marginal tax rate or a change in skill, the maximum remains global.

Along the nonlinear income tax schedule, we define the compensated elasticity of earnings with respect to the marginal retention rate 1 -T (.) as the elasticity of earnings for individuals of type (w, θ) to a change in the marginal retention rate by a constant amount τ, while leaving unchanged the level of tax at y = Y(w, θ), i.e.:

ε(w, θ) def ≡ 1 -T (Y(w, θ)) Y(w, θ) ∂Y ∂τ c (19a)
which is positive as shown in Appendix A.4 and where the superscript "c" stands for "compensated".

Along the nonlinear income tax schedule, the income effect is defined as the behavioral response to a lump-sum change m in tax liability:

η(w, θ) def ≡ ∂Y ∂m i (19b)
where the superscript "i" stands for "income effect". We have η(w, θ) < 0 if leisure is a normal good (see appendix).

Let us use the dot above a variable for the partial derivative of this variable with respect to skill w. One can define the elasticity α(w; θ) of earnings with respect to the skill level:

α(w, θ) def ≡ w Y(w, θ) Ẏ(w, θ). ( 19c
)
which is positive thanks to Assumption 2. We can note that ε(w, θ), η(w, θ) and α(w, θ) denote total responses of earnings since they take into account the nonlinearity of the tax schedule as in [START_REF] Jacquet | Optimal redistributive taxation with both extensive and intensive responses[END_REF], see also Scheuer and Werning (2016).

We define the social marginal welfare weights associated with workers of type (w, θ) expressed in terms of public funds by:

g (w, θ) def ≡ Φ U (U (w, θ) ; w, θ) U c (C(w, θ), Y(w, θ); w, θ) λ . ( 20
)
The government values giving one extra dollar to a worker (w, θ) as a gain of g(w, θ) dollar(s) of public funds.

Let h(y|θ) denote the conditional income density within group θ at income y and H(y|θ) def ≡ y 0 h(z|θ)dz the corresponding conditional income CDF. According to (19c) and Assumption 2, income Y(•, θ) is strictly increasing in skill within each group. We then have H (Y(w, θ)|θ) ≡ F(w|θ) for each skill level w. Differentiating both sides of this equality with respect to w and using (19c), the two densities are linked by:

h (Y(w, θ)|θ) = f (w|θ) Ẏ(w, θ) ⇔ Y(w, θ) h (Y(w, θ)|θ) = w f (w|θ) α(w, θ) . ( 21
)
The unconditional income density at income Y(w, θ 0 ) is given by: ĥ(Y(w, θ 0 ))

def ≡ θ∈Θ h(Y(W(w, θ), θ)|θ) dµ(θ). (22a)
The mean (total) compensated elasticity at income level y = Y(w, θ 0 ) is:

ε(Y(w, θ 0 )) = θ∈Θ ε (W(w, θ), θ) h(Y(W(w, θ), θ)|θ) ĥ(Y(w, θ 0 )) dµ(θ). (22b)
where each within-group compensated elasticity of earnings is timed by the relative proportion h(y|θ)/ ĥ(y) of individuals in the corresponding group among individuals who earn y. The mean (total) income effect at income level Y(w, θ 0 ) is:

η(Y(w, θ 0 )) = θ∈Θ η (W(w, θ), θ) h(Y(W(w, θ), θ)|θ) ĥ(Y(w, θ 0 )) dµ(θ). (22c)
Finally, the mean marginal social welfare weight at income level y = Y(w, θ 0 ) is:

ĝ(Y(w, θ 0 )) = θ∈Θ g (W(w, θ), θ) h(Y(W(w, θ), θ)|θ) ĥ(Y(w, θ 0 )) dµ(θ). (22d)
As shown in Appendix A.4, we can now rearrange the first-order conditions (14a) and (14b) displayed in Proposition 1 to obtain the optimal marginal tax rate in terms of sufficient statistics, more precisely in terms of the mean compensated elasticity, mean income effect, mean marginal social weights and the unconditional income density.

Proposition 2. Under assumption 2, with a twice differentiable T(.), with the second order condition of the individual maximization program that holds strictly and the function y → U (y -T(y), y; w, θ)

that admits a unique global maximum, the optimal tax schedule satisfies:

T (y) 1 -T (y) = 1 ε(y) • 1 -Ĥ(y) y ĥ(y) • 1 - ∞ y [ ĝ(z) + η(z) • T (z)] • ĥ(z)dz 1 -Ĥ(y) (23a) 1 = ∞ 0 ĝ(z) + η(z) • T (z) • ĥ(z)dz. (23b)
The optimal tax rate given in Equation (23a) depends now on three terms that encapsulate sufficient statistics only: (a) the behavioral responses to taxes denoted by 1/ε(y), which, in the vein of [START_REF] Ramsey | A Contribution to the Theory of Taxation[END_REF], is the inverse of the mean compensated elasticity; (b) the social preferences and income effects 1 -∞ y ĝ(z) + η(z) T (z) ĥ(z)dz / (1 -Ĥ(y)), which indicates the distributional benefits of increasing the tax liability by one unit for all workers with incomes above y and (c) the shape of the income distribution measured by the inverse of the local Pareto parameter (1 -Ĥ(y))/(y ĥ(y)) of the income distribution. In Equation (23a), the sufficient statistics' optimal tax formula of [START_REF] Saez | Using Elasticities to Derive Optimal Income Tax Rates[END_REF] is generalized to the multidimensional context.

Equations ( 23a) and (23b) can also be obtained, in a straightfoward way, using the tax perturbation approach as in [START_REF] Jacquet | Optimal redistributive taxation with both extensive and intensive responses[END_REF]. Thanks to these equations, which characterize the optimal tax schedules in terms of incomes, we know the weight that multiplies each sufficient statistic. In (23a), the unconditional income density ĥ(Y(w, θ 0 )) is the weight that multiplies each mean marginal social welfare weight and each mean income effect. In contrast, the weights that are used in the optimal tax formulas in terms of skills (i.e. Equations (14a) and (14b)) are distinct. These weights are the conditional skill densities times the corresponding skill. 7

7 More precisely, in the left-hand side of Equation (14a), the term - [w,θ] which is equal to the ratio of ε(w, θ) and α(w, θ) (see Equation (35) in the appendix), is weighted by the conditional density times the skill, W(w, θ) f (W(y, θ)|θ). And, in the right-hand side of (14a), which encapsulates the mechanical and income effects, the weights are the conditional skill densities.

v y [w,θ] w•v yw

V Signing optimal marginal tax rates

In the literature, it has been highlighted that, with multidimensional heterogeneity, negative marginal tax rates can be optimal. In [START_REF] Cuff | Optimality of workfare with heterogeneous preferences[END_REF], Boadway et al. (2002a), [START_REF] Brett | How optimal nonlinear income taxes change when the distribution of the population changes[END_REF], [START_REF] Choné | Negative Marginal Tax Rates and Heterogeneity[END_REF] and [START_REF] Lockwood | De Gustibus non est Taxandum: Heterogeneity in Preferences and Optimal Redistribution[END_REF], individuals differ along their skills and preferences for leisure, and the social planner has weighted utilitarian preferences. In this context, individuals who pool at the same income level are characterized by different marginal social welfare weights, Φ U (U(w, θ); w, θ)u (C(w, θ))/λ. Therefore, the mean marginal social welfare weight, θ∈Θ,w∈R + Φ U (U(w,θ);w,θ)u (C(w,θ)) λ f (w|θ)dw dµ(θ), may not be decreasing with income which opens the way to negative marginal tax rates. One may wonder whether heterogeneous behavioral responses for individuals who earn the same level of income could also open the way to negative marginal tax rates. We now show that this is not the case (if the population earning a given income level remains homogeneous in terms of social welfare weights).

Proposition 3. Under utilitarian or maximin social preferences, optimal marginal tax rates are positive.

Proposition 3 generalizes to the multidimensional case [START_REF] Mirrlees | An Exploration in the Theory of Optimum Income Taxation[END_REF]'s result of positive optimal tax rates (which was obtained under additively separable preferences). 8 The proof, which relies on Proposition 1, can be found in Appendix A.5. Proposition 3 emphasizes that optimal marginal tax rates are positive as soon as all individuals who earn the same income y are characterized (i) by the same marginal utility of consumption, which is ensured by the additive separability assumption and (ii) by the same marginal social welfare Φ u , which is ensured by utilitarian or maximin social objective. In such a case, all individuals who earn the same income are characterized by the same welfare weights. Therefore, the cause of negative marginal tax rates due to decreasing mean marginal social welfare weights emphasized in [START_REF] Cuff | Optimality of workfare with heterogeneous preferences[END_REF] and follow-up papers does not apply.9

VI Concluding Comments

In this paper, we have proposed an allocation perturbation method to derive optimal nonlinear income tax schedules when taxpayers differ along several characteristics and when taxation cannot be conditioned on them. In this context, we have shown that marginal tax rates are positive under utilitarian preferences and maximin.

To illustrate the generality of our results in this concluding section, we now provide alternative tax problems that one can easily solve in our framework. For each of them, we explain what y, w, θ represent so that the interpretation of the results is straightforward.

Optimal joint taxation of labor and non-labor income

Consider individuals that have two sources of taxable income: a non-labor income z and a labor income yz. Those incomes are jointly taxed and the tax function does not distinguish between both incomes. This applies, for instance, in countries like France where income received from renting property and entrepreneurial income are jointly taxed with labor income.

As explained in [START_REF] Scheuer | Adverse selection in credit markets and regressive profit taxation[END_REF], a single nonlinear tax schedule is also the system that is in place for employed workers and self-employed small business owners in many countries, including the U.S.. In this case, y is the total taxable income and we interpret θ as the ability to earn non-labor income z and w as the skill. Individuals of type (w, θ) solve: max y,z U (y -T(y), yz, z; w, θ)

where two decision variables appear instead of one variable in the core of our paper. This program can be solved sequentially, the first step being the choice of non-labor income z for a given taxable income y which leads to U (c, y; w, θ) def ≡ max z U (c, yz, z; w, θ). The second step is the choice of y as in Equation ( 1). In the process, one simply needs to ensure the semiindirect utility function U (•, •; w, θ) verifies Assumption 2.

Optimal joint income taxation of couples

The joint income taxation of couples is a variant of the previous application, in which yz is the labor income of one individual and z is the one of his/her partner. The tax does not distinguish between yz and z and only depends on the sum of both incomes, y (as in France, Germany and the US). We redefine w and θ as the respective skill level of each member of the couple. The optimal tax schedules derived in this paper are then interpreted as the optimal tax schedules when the couple is the tax unit and each partner decides along the intensive margin.

So far, previous attempts in the literature [START_REF] Kleven | The Optimal Income Taxation of Couples[END_REF][START_REF] Cremer | Income taxation of couples and the tax unit choice[END_REF]) have stopped short of obtaining these nonlinear tax schedules.

Optimal income taxation with tax avoidance

In this application, w is the skill and θ is the ability to avoid taxation. We assume that tax enforcement (penalty, monitoring, etc.) is given. We denote z the sheltered labor income (i.e. income that is not taxed at all) and y + z the (total) labor income. The tax only depends on the taxable income y. Consumption becomes c + z, with c = y -T(y) being the after-tax income.

All results obtained in this paper are valid in this context when one simply makes sure that Assumption 2 holds.

A Appendix

A.1 Proof of Lemma 4

Proof The proof consists of two steps. In step (i), we show that there exists at most one incentive-compatible allocation (w, θ) → (C(w, θ), Y(w, θ)) that verifies Assumption 3 and such that (C(w, θ 0 ), Y(w, θ 0 )) = (C(w, θ 0 ), Y(w, θ 0 )). In step (ii), we show that this allocation verifies the whole set of incentive constraints (6).

Step (i). To build up the entire incentive-compatible allocation (w, θ) → (C(w, θ), Y(w, θ)), we must choose (C(w, θ 0 ), Y(w, θ 0 )) = (C(w, θ 0 ), Y(w, θ 0 )) at any skill level. For each group θ, Y(•, θ) verifies Assumption 3 if and only if its reciprocal Y -1 (•; θ) is smoothly increasing. Let y ∈ R + be an income level. As Y(•, θ 0 ) is smoothly increasing from Assumption 3, there exists a unique skill level w such that y = Y(w, θ 0 ). Then according to Lemma 3, among individuals of group θ, only those of skill W(w, θ) must be assigned to the income level y = Y(w, θ 0 ) to verify incentive-compatibility.10 Therefore, Y -1 (•, θ) must be defined by:

Y -1 (•, θ) : y Y -1 (•,θ 0 ) -→ w = Y -1 (y, θ 0 ) W(•,θ) -→ Y -1 (y, θ).
Y -1 (•, θ) is then smoothly increasing as a combination of two smoothly increasing functions. Moreover, since for each type (ω, θ), there exists a single skill level ω such that Y(ω, θ) = Y(w, θ 0 ), incentive compatibility requires that C(ω, θ) also needs to be equal to C(w, θ 0 ). This ends the proof of step (i).

Step (ii). Note that the allocation (w, θ) → (Y(w, θ), C(w, θ)) is built in such a way that one has Y(ω, θ) = Y(w, θ 0 ) and C(ω, θ) = C(w, θ 0 ) if and only if ω = W(w, θ) and ( 10) holds. Differentiating in w both sides of these two equations and rearranging terms, we obtain

Ċ (w, θ 0 ) Ẏ (w, θ 0 ) = Ċ (W(w, θ), θ 0 ) Ẏ (W(w, θ), θ 0 ) . As w → (C(w, θ 0 ), Y(w, θ 0 )
) is assumed to verify the within-group incentive constraints in Equation (8b), we know that the left-hand side of the above equation is equal to M (C(w, θ 0 ), Y(w, θ 0 ); w, θ 0 ).

Using the definition of W(•, θ), we have that w → (C(w, θ), Y(w, θ)) also verifies Equation (8b). From Lemma 2, it thus verifies the within-group incentive constraints (7). We now check whether the inequality ( 6) is verified for any (w, w , θ, θ ) ∈ R 2 + × Θ 2 . We know there exists ω ∈ R + such that Y(ω, θ) = Y(w , θ ) and C(ω, θ) = C(w , θ ). The incentive constraints in (6) are therefore equivalent to:

U (C(w, θ), Y(w, θ); w, θ) ≥ U (C(ω, θ), Y(ω, θ); w, θ) .
The latter inequality is verified as w → (C(w, θ), Y(w, θ)) satisfies Equation (8b).

A.2 Derivation of Equation (17)

Proof

To derive (17), we must compute the various Gâteaux derivatives at t = 0. For A = C, Y, U and a given δ, the Gâteaux derivative of A in the direction ∆ Y (•, •; δ) at t = 0 is denoted Â(x, θ; δ). Let us remind its definition: Â(x, θ; δ) υ yw (Y(ω, θ); ω, θ) Ŷ(ω, θ; δ) dω. (24d)

Using Equations ( 12), (24a) and (24c), the Gâteaux derivative of the Lagrangian ( 16) is: υ yw (Y(x, θ); x, θ) Ŷ(x, θ; δ) dx × ∞ W(w,θ)

∂ L ∂t (0; δ) =
Φ U [x, θ] λ - 1 u [x, θ] f (x|θ)dx dµ(θ).
Dividing the first-order condition ∂ L ∂t (0; δ) = 0 by w w-δ υ yw (Y(x, θ 0 ); x, θ 0 ) Ŷ(x, θ 0 ; δ) dx im- plies, using (24b) and (24d), that

Φ U [x, θ] λ - 1 u [x, θ]
x W(w-δ,θ) υ yw (Y(x, θ); x, θ) Ŷ(x, θ; δ) dx W(w,θ) W(w-δ,θ) υ yw (Y(x, θ); x, θ) Ŷ(x, θ; δ) dx f (x|θ)dx + ∞ W(w,θ)

Φ U [x, θ] λ - 1 u [x, θ] f (x|θ)dx    dµ(θ).
We finally take the limit of the latter equality when δ tends to 0. Let us consider the first term in the right-hand side of (26). Since

A.3 Proof of Lemma 5

With one-dimensional heterogeneity, we only consider within-group incentive constraints. Adopting a first-order approach, only (8a) is considered when building up the Hamiltonian: Y(w, θ) -C (Y(w, θ), U(w, θ); w, θ) + Φ (U(w, θ); w, θ) λ • f (w|θ)q(w|θ) • v w (Y(w, θ); w, θ) .

where Y(w, θ) and U(w, θ) are the control and state variables respectively. Using (12), the necessary conditions are:

0 = 1 - v y [w, θ] u [w, θ] • f (w|θ) -q(w|θ) • v yw [w, θ] (28a) -q (w|θ) = Φ U [w, θ] λ - 1 u [w, θ]
• f (w|θ) (28b) 0 = q(0|θ) (28c) 0 = lim w →∞ q(w|θ).

Combining ( 28b) with (28d) leads to

q(w|θ) = ∞ w Φ U [w, θ] λ - 1 u [w, θ]
• f (ω|θ)dω.

(28e)

Combining ( 3), ( 2), (28a) and (28e) leads to (18a). Combining (28c) with (28e) leads to (18b).

A.4 Proof of Proposition 2

Define a reform of a tax schedule y → T(y) with its direction, which is a differentiable function y → R(y) defined on R + , and with its algebraic magnitude m ∈ R. After a reform, the tax schedule becomes y → T(y)m R(y) and the utility of an individuals of type (w, θ) is:

  θ; t, δ) -A(w, θ) t .By definition we get: Ŷ(x, θ 0 ; δ) = ∆ Y (x; δ), and from (15b) we obtain:Ŷ(x, θ; δ) = 0 if x ∈ [0, W(wδ, θ)] ∪ [W(w, θ), +∞) . (24a)Equations (15c) imply that the Gâteaux derivatives of utilities are nil for skill below W(wδ, θ).For skills x between W(wδ, θ) and W(w, θ), Equation (15e) implies:Û(x, θ; δ) = -x W(w-δ,θ) υ yw (Y(ω, θ 0 ); ω, θ 0 ) Ŷ(ω, θ 0 ; δ) dω. (24b)For skill x above W(w, θ), according to (15f), we have:Û(x, θ; δ) = -W(w,θ) W(w-δ,θ) υ yw (Y(ω, θ 0 ); ω, θ 0 ) Ŷ(ω, θ 0 ; δ) dω.(24c)Moreover, Equation (15h) implies that the Gâteaux derivatives of income must verify: w w-δ υ yw (Y(ω, θ 0 ); ω, θ) Ŷ(ω, θ 0 ; δ) dω =

  y (Y(x, θ); x, θ) u (C(x, θ)) Ŷ(x, θ; δ) f (x|θ)dx (25) + W(w,θ) W(w-δ,θ) Φ U [x, θ] λ -1 u [x, θ] Û(x, θ; δ) f (x|θ)dx -

  δ,θ) 1 -υ y (Y(x, θ); x, θ) u (C(x, θ)) Ŷ(x, θ; δ) f (x|θ)dx W(w,θ) W(w-δ,θ) υ yw (Y(x, θ); x, θ) Ŷ(x, θ; δ)

  δ,θ) υ yw (Y(x, θ); x, θ) Ŷ(x, θ; δ) dx W(w,θ) W(w-δ,θ) υ yw (Y(x, θ); x, θ) Ŷ(x, θ; δ) dx ∈ [0, 1] we get that: θ∈Θ W(w,θ) W(w-δ,θ) Φ U [x, θ] λ -1 u [x, θ]

  x W(w-δ,θ) υ yw (Y(x, θ); x, θ) Ŷ(x, θ; δ) dx W(w,θ) W(w-δ,θ) υ yw (Y(x, θ); x, θ) Ŷ(x, θ; δ) dx f (x|θ)dxdµ(θ) ≤ θ∈Θ W(w,θ) W(w-δ,θ) Φ U [x, θ] λ -1 u [x, θ]f (x|θ)dxdµ(θ) .

  -T(y) + m R(y))υ(y; w, θ).(29) We denote by Y R (m; w, θ) the income of workers of types (w, θ) after the reform and her con-sumption becomes C R (m; w, θ) = Y R (m; w, θ) -T(Y R (m; w, θ)) + m R(Y R (m; w, θ)). When m = 0, we have Y R (0; w, θ) = Y(w, θ) and C R (0; w, θ) = C(w, θ).Applying the envelope theorem to (29), we get:∂U R ∂m (m; w, θ) = u c C R (m; w, θ) R(y).(30)Using (3), the first-order condition associated to (29) equalizes to zero the following expression:Y R (y, m; w, θ) def ≡ 1 -T (y) + m R (y) -M (y -T(y) + m R(y), y; w, θ) . (31)For simplicity, we drop the superscript R and write Y y (Y(w, θ); w, θ) for Y R y (Y(w, θ), 0; w, θ). 11 We define behavioral responses to tax reforms of direction R by applying the implicit function theorem to Y R (y, m; w, θ) = 0 at m = 0, which yields: ∂Y R ∂m (0; w, θ) = -Y R m (Y(w, θ), 0; w, θ) Y y (Y(w, θ), 0; w, θ) (32)

John Weymark has largely contributed to the literature on optimal income taxation, see e.g.[START_REF] Weymark | A reduced-form optimal nonlinear income tax problem[END_REF] andWeymark (1987) and[START_REF] Brett | How optimal nonlinear income taxes change when the distribution of the population changes[END_REF].

Our paper studies the optimal tax system when individual characteristics, despite being observable by the tax authority, cannot be used as tags[START_REF] Akerlof | The Economics of "Tagging" as Applied to the Optimal Income Tax, Welfare Programs, and Manpower Planning[END_REF], due to legal and/or horizontal equity reasons.

Our definition of "group" is identical to the one inWerning (2007, p.13).

A smoothly increasing (decreasing) function is also called an increasing (decreasing) diffeomorphism for which the derivative maps the positive real line onto itself.

InJacquet and Lehmann (2020, Proposition 5), we show that the assumption of a smoothly-increasing-in-types allocation amounts to assuming: (i) twice differentiability of the tax function T(•), that (ii) for all (w, θ) ∈ R * + × Θ, the second-order condition associated to the individual maximization program holds strictly and that (iii) for all (w, θ) ∈ R * + × Θ, the function y → U (y -T(y), y; w, θ) admits a unique global maximum over R + .

For instance, we never found cases where the second-order incentive-compatibility constraints were violated in the large set of simulations we run on US data with taxpayers differing in terms of gender and labor supply elasticities, see[START_REF] Jacquet | Optimal redistributive taxation with both extensive and intensive responses[END_REF].

[START_REF] Hellwig | A contribution to the theory of optimal utilitarian income taxation[END_REF], under a utilitarian criterion, positive optimal tax rates are obtained with more general preferences.

 9 If the utility function u(•) in (1) were parameterized by type w and θ while v(•) were simply parameterized by w, individuals who earn the same income would have distinct social marginal welfare weights. This could drive negative marginal tax rates. Similarly, if both u(•) and v(•) were parameterized by w and θ, one would also expect negative marginal tax rates. Let us stress that our method could not be used in this framework since the pooling function (10) cannot depend simultaneously on Y and C.

Hence function W(•, θ) coincides with the pooling function W(•, θ).

Indeed, at m = 0, Y R y does no longer depend on the direction R of the tax reform.

She gratefully acknowledges support by Labex MME-DII. The usual disclaimer applies.

As the right hand-side of the latter inequality tends to 0 when δ tends to 0, the limit of (26) when δ tends to zero leads to:

W(w-δ,θ) υ yw (Y(x, θ); x, θ) Ŷ(x, θ; δ) dx dµ(θ)

= θ∈Θ,x≥W(w,θ)

By continuity, the variations of f (x|θ), υ y (Y(x, θ); x, θ), υ yw (Y(x, θ); x, θ) and u (c(x, θ)) within the skill intervals [W(wδ, θ), W(w, θ)] are of second-order when δ tends to 0. As Θ and intervals [W(wδ, θ), W(w, θ)] are compact, for any small e > 0, there always exists δ(e) such that for all (x, θ) ∈ [W(w -δ(e), θ), W(w, θ)] × Θ, one has:

and

and therefore, for all δ < δ(e):

Hence, left-hand side of ( 27) is equal to the left-hand side of ( 17).

where:

Using (2) and plugging R(Y(w, θ)) = 0 and R (Y(w, θ)) = 0 into (33b), the compensated elasticity of earnings (19a) can be rewritten as:

, the income effect (19b) can be rewritten as:

which is negative if leisure is a normal good, since then M c > 0. The elasticity α(w; θ) of earnings with respect to the skill level can be expressed as:

Dividing (34a) by (34c) we get:

Plugging (34a) into (34b) leads to:

It is then straightforward to obtain:

Let y ∈ R + . Since Y y (Y(w, θ); w, θ) < 0, there exists a single skill level w such that y = Y(w, θ 0 ). From (2), we know that:

The term in the left-hand side integral of (14a) can be rewritten as:

The first equality is obtained using Equation ( 35). The second equality uses (21). It implies with (22b) that Equation (14a) can be rewritten as:

where J(w) is defined by the right-hand side of (14a). J(•) admits for derivative J(w) where:

where ( 20) has been used. Deriving with respect to the skill w both sides of (9) and of C(w, θ 0 ) = Y(w, θ 0 ) -T (Y(w, θ 0 )), we obtain:

We thus obtain:

Using ( 21) and (38), J(w) can be rewritten as:

Using (36) and (22d), we get:

As J(w) = x≥w (-J(x))dx, we get

Changing variables by posing z = Y(x, θ 0 ), we get

Plugging ( 39) into ( 38) gives (23a). Combining (14b) and (39) leads to (23b).

A.5 Proof of Proposition 3

Let us denote

the ratio of the right-hand side of (14a) at the skill level w divided by u (Y(w, θ 0 ) -T(Y(w, θ 0 ))).

According to Proposition 1, Equation (14a) and υ y > 0 > υ yw , the sign of T (Y(w, θ 0 )) is the sign of K(w).

Under utilitarian preferences, Φ u = 1. Changing variable in (40) from x to t such that x = W(t, θ), (i.e. Y(x, θ) ≡ Y(t, θ 0 ) and C(x, θ) ≡ C(t, θ 0 ) according to ( 9)), we get:

The derivative of K(w) has the sign of 1/λ -1/u (C(w, θ 0 )), which is decreasing in w because of the concavity of u(•). Moreover, lim w →∞ K(w) = 0 and Equation (14b) implies that K(0) = 0.

Therefore, K(•) first increases and then decreases. It is thus positive for all (interior) skill levels. So, optimal marginal tax rates are positive. Under maximin, one has U(x, θ) > U(0, θ) for all x > 0 from (8a). Therefore, within each group, the most deserving individuals are those whose skill w = 0. The maximin objective implies Φ U [x, θ] = 0 for all x > 0. Hence, Equation (40) simplifies to:

for all w > 0, which is always positive, thereby leading to positive optimal marginal tax rates.