Liang Z Tan

Milan Orlita

Marek Potemski

James Palmer

Claire Berger

Walter A De Heer

Steven G Louie

Gérard Martinez

SU(4) Symmetry Breaking Revealed by Magneto-optical Spectroscopy in Epitaxial Graphene

Refined infra-red magneto-transmission experiments have been performed in magnetic fields B up to 35 T on a series of multi-layer epitaxial graphene samples. Following the main optical transition involving the n = 0 Landau level, we observe a new absorption transition increasing in intensity with magnetic fields B ≥ 26T. Our analysis shows that this is a signature of the breaking of the SU(4) symmetry of the n = 0 LL. Using a quantitative model, we show that the only symmetry breaking scheme consistent with our experiments is a charge density wave (CDW).

I. INTRODUCTION

In multicomponent quantum Hall systems, interaction effects lead to a rich variety of broken symmetry ground states. In graphene, the spin and valley degrees of freedom of the lowest Landau level (LL) form an SU(4) symmetric quartet. Refined transport experiments have shown evidence of a broken symmetry state [1][2][3][4][5][6][7][8], but there is no clear consensus on its nature [9][10][11][12], and how it is affected by different substrates and disorder. While the spin degree of freedom has been probed in tilted magnetic fields [1,2,5,6], we show here that the valley degree of freedom can be accessed by examining the signatures of optical phonons in magneto-transmission spectra. In this paper, we show that our observation of a new absorption transition supports the existence of a charge density wave (CDW) in our epitaxial graphene samples.

In a quasiparticle picture, charge carriers in graphene are characterized by a Dirac-like spectrum around the K and K ′ equivalent points ("valley") of the Brillouin zone of the hexagonal crystal lattice. As a consequence, the application of a magnetic field B perpendicular to the plane of the structure splits the electronic levels into Landau levels (LL) indexed by n, with specific energies E n = sgn(n)v F 2e B|n| where n are integers including 0 (v F being the Fermi velocity). In this paper, we are concerned with how a broken symmetry phase can be observed in infra-red magneto-optical transitions involving the n = 0 LL (i.e., transitions from n = -1 to n = 0 or from n = 0 to n = 1 equivalent to a cyclotron resonance (CR) transition in the quantum limit) with an energy E 01 = v F √ 2e B [13]. Our previous work [14] has reported on the magnetic field dependence of this transition revealing its interaction with the K-phonon. Besides this specific interaction, we observed that the basic broadening γ 01 (B) ∝ √ B of the transition had an ad-ditional component proportional to B in contrast to all theoretical models [15]. This could be already a sign of the breaking of the valley degeneracy.

In the present work, we use the Γ-phonon at the Brillouin zone center as a probe of the valley symmetry breaking. In the absence of valley symmetry breaking, the Γ-phonon does not affect the infra-red absorption spectrum because the electron-phonon matrix elements are of opposite signs for the K and K ′ valleys [16]. However, one expects to see signs of valley symmetry breaking when the energy E 01 (B) is larger than that of the optical Γ-phonon ( ω Γ = 0.196eV). It turns out, indeed, that when that condition is reached, a new optical transition develops at an energy higher than the main line (Fig. 1). We interpret this as a signature of the breaking of the SU(4) symmetry. A model has been established to reproduce these findings and applied to the different phases which have been proposed.

In Sec. II, we discuss these experimental observations and methods in more detail. We first interpret our experimental findings within a simplified model for valley symmetry breaking in Sec. III before deriving a more complete Hamiltonian in Sec. IV and calculating the optical conductivity in various broken symmetry phases in Sec. V. Finally, a comparison between experiment and theory is presented in Sec. VI, followed by conclusions in Sec. VII.

II. EXPERIMENTAL OBSERVATIONS AND METHODS

In our experiment, precise infra-red transmission measurements were performed on multi-layer epitaxial graphene samples, at 1.8 K, under magnetic fields up to 35 T. The light (provided and analyzed by a Fourier transform spectrometer) was delivered to the sample by means of light-pipe optics. All experiments were performed with nonpolarized light, in the Faraday geometry with the wave vector of the incoming light parallel to the magnetic field direction and perpendicular to the plane of the samples. A Si bolometer was placed directly beneath the sample to detect the transmitted radiation. The response of this bolometer is strongly dependent on the magnetic field. Therefore, in order to measure the absolute transmission TA (B, ω), we used a sample-rotating holder and measure for each value of B a reference spectrum through a hole. These spectra are normalized in turn with respect to TA(0, ω) to obtain a relative transmission spectrum TR(B, ω) which only displays the magnetic field dependent features. Those spectra are presented in Fig. 2.

The samples were grown [17] on the C-terminated surface of SiC and display the characteristic transmission spectra of isolated graphene monolayers that arise from rotational stacking of the sheets [18]. The thickness d of the SiC substrate has to be reduced significantly in order to minimize the very strong double-phonon absorption of SiC in the energy range of interest. In the first series d was reduced to 60µm and related samples have been used to perform the experiments reported earlier [14]. One of them, named S4, was used to compare the data with those obtained on sample S5 from a new series where the thickness d was further reduced down to 32µm. We compare in Fig. 2 the transmission spectra, at high fields, for samples S4 and S5. Technically speaking, the optical response of both samples is almost the same, showing that they have a similar number of active layers.

Taking into account all layer dielectric properties of each sample in a multi-layer dielectric model, we determine the effective number, N ef f , of graphene sheets with their respective carrier density (see Appendix). For samples S4 and S5, we have found that N ef f = 7 with carrier densities {5.5, 3.2, 1.8, 1.1, 0.5, 0.1, 0.1}10 12 cm -2 and {6.0, 3.5, 2.0, 1.2, 0.5, 0.1, 0.1}10 12 cm -2 respectively. These carrier densities are fixed for each sample.

The transmission spectra of sample S5 at high magnetic fields are displayed in Fig. 1. We observe a new transition occurring at an energy higher than that of the CR line, growing in intensity when increasing the magnetic field. This behavior cannot be explained without breaking the SU(4) symmetry in graphene. In order to characterize more clearly these findings, one can treat the data, as a first step and in a very rough way, extracting from the transmission data the real part of the effective diagonal component of the conductivity σ xx (ω, B) [20]. We have deconvoluted this result with two Lorentzians of equal width, extracting the evolution of the two extrema with the magnetic field. The resulting energies are displayed in Fig. 3 (top panel) for samples S4 and S5.

Though the procedure adopted at this initial level is quite rough, it provides important information: (i) The evolution of the lower energy line varies at low fields like B 1/2 (function F2(B) in Fig. 3) with a coefficient proportional to the Fermi velocity v F and ends at higher fields with a similar dependence (function F1(B)) but with a smaller value of v F which is, by itself, a sign of some interaction occurring at an energy close to that of the Γ-phonon; (ii) The second component of the deconvolution always appears at energies larger than that of the Γ-phonon; (iii) In principle, in the SU(4) symmetric picture, it is not possible to explain the occurrence of an additional transition, growing in intensity with B, at higher energies than the main transition line; (iv) It is therefore clear that the Γ-phonon plays a crucial role though it should not, indicative that the SU(4) symmetry is broken. Using these observations, we now have some guidelines to develop a theory which can explain quantitatively the experimental observations. In addition, we note that results for samples S4 and S5 are quite similar within the experimental errors. Knowing that the active layers which contribute to the E 01 transition should have a filling factor ν ≤ 2 ( ν = N s Φ 0 /B, Φ 0 being the flux quantum and N s the carrier density) and that, in samples S4 and S5, the carrier density for active layers do not have the same sequence, the physical mechanisms describing the experimental findings should not be very dependent on the doping of active layers. This is indeed the case as discussed below.

III. SIMPLIFIED MODEL FOR VALLEY SYMMETRY BREAKING

To illustrate how the electron-phonon interaction and the valley symmetry breaking give rise to the observed features in the transmission spectrum, we first introduce a simplified model for the interaction of the Γ phonon with the E 01 excitation, before discussing the full SU(4) calculation. The simplified model provides a minimal description of the valley symmetry breaking by neglecting the spin degree of freedom in the n = 0 LL. We assume that K and K ′ sublevels of the n = 0 LL are separated in energy by ∆ V , and have different filling factors ν K and ν K ′ . Considering just the n = 0 to n = 1 transitions, the interaction with the Γ phonon is captured by the Hamiltonian (in the basis of creating an electronic excitation in K, electronic excitation in K ′ , and a Γ phonon, in that order)

H =   E 01 -∆ V /2 0 g ph √ ν K 0 E 01 + ∆ V /2 -g ph √ ν K ′ g ph √ ν K -g ph √ ν K ′ ω Γ   . ( 1
)
where g ph characterizes the electron-Γ-phonon interaction. The optical conductivity is calculated using the Green's function formalism introduced by Toyozawa [21].

The diagonal component of the conductivity is:

σ xx ( ω) = 1 ω ImM † x G( ω)M x (2)
where the Green's function is G = ( ω -H -iη) -1 , with η → 0 + (see Sec. IV B). The optical matrix elements for the simplified model are

M † x = √ ν K , √ ν K ′ , 0 .
The simplified model explains the splitting of the main transition line in the two limits E 01 ≈ ω Γ and E 01 ≫ ω Γ . In the absence of valley splitting (∆ V = 0 and ν K = ν K ′ ), the eigenstates of the pure electronic part of H in Eq. 1 are valley-symmetric and valley-antisymmetric combinations of E 01 transitions, i.e.

1 2 c † 1,K c 0,K + c † 1,K ′ c 0,K ′ |GS and 1 2 c † 1,K c 0,K -c † 1,K ′ c 0,K ′ |GS , respectively
, where |GS denotes the ground state and c † n,K are creation operators at LL n and valley K. The valley-symmetric combination is infra-red active but does not interact with the Γ phonon, while the valley-antisymmetric combination is infra-red inactive and interacts with the Γ phonon. The symmetry breaking valley splitting term ∆ V allows both eigenmodes to interact with the Γ phonon while remaining infra-red active, inducing a splitting of the main transmission line in the vicinity of the Γ phonon frequency. Away from the Γ phonon frequency (E 01 ≫ ω Γ ), E 01 transitions at K and K ′ interact weakly with the phonon and the splitting of the main transmission line is controlled directly by the energy difference ∆ V .

IV. THEORY OF MAGNETO-PHONON RESONANCE IN THE PRESENCE OF SU(4) SYMMETRY BREAKING

We reintroduce the spin degree of freedom and the n = -1 to n = 0 transitions in order to obtain a quantitative understanding of the experiment. We consider different theoretical models of the n = 0 LL SU(4) symmetry breaking, taking into account the effects of ν = 0 and disorder by introducing Gaussian broadening into a mean field theory (Sec. IV A). Different symmetry-breaking phases are represented in the mean field theory by different orderings and filling factors of the four sublevels of the n = 0 LL. We consider four candidate symmetrybreaking phases that have been proposed in the literature [9]: Ferromagnetic(F), Charge Density Wave (CDW), Canted Antiferromagnetic (CAF) and Kekulé-distortion (KD), and calculate the optical conductivity using Eq. 2 with the appropriate Hamiltonian H for each phase. Treating these phases on the same footing (detailed in Sec. V), we find that each phase results in characteristic features in the evolution of the transmission spectrum as a function of the magnetic field. By examining the intensities and positions of the transmission lines, we identify the symmetry broken phase in the samples used in our experiment as the CDW type [27,28].

A. Description of the ground state

We assume that the ground state is a single Slater determinant of the form:

|GS = 4 j=1 Nj mj =1 Ψ † j,mj |0 (3)
the index j runs over the 4-dimensional spin/valley space and m j describe the "guiding center" degree of freedom. The state (j, m j ) is represented by the wavefunction ξ j φ mj ( -→ r ) where ξ j is a four-component spinor and φ mj ( -→ r ) is the orbital part of the wavefunction. These wavefunctions belong to the n = 0 landau level (LL) of graphene. The occupation numbers N j count the number of j states that are occupied in this ground state.

There are different models proposed to describe the symmetry-broken phase of graphene which have been reviewed by Kharitonov [9]. For a given model, we assume that the system is polarized along a certain direction in j-space. For instance, with increasing order of energies, j = 1, 2, 3, 4 corresponds to (K ′ ↑, K ′ ↓, K ↑, K ↓) in the charge density wave (CDW) phase. The remaining degrees of freedom, φ mj ( -→ r ) and N j , are treated as variational parameters, subject to the constraint

N 1 + N 2 + N 3 + N 4 = N .
We minimize the energy of the ground state E GS = GS|H 0 + H e-e + H disorder |GS . Here H 0 is the single part of the Hamiltonian without disorder, H e-e the interaction term and H disorder the disorder potential. Because we assume a single Slater determinant, we can apply mean-field theory and obtain single-particle energy levels E j,mj (The origin of the energies is taken to be at the energy of the n = 0 LL of the non-interacting system).

In a system with finite disorder, the energy levels E j,mj are clustered about mean values E j = avg mj E j,mj . We remove the m j degrees of freedom by replacing the energy levels E j,mj by broadened energy levels centered at E j .

There is a Fermi level E F which fixes the occupation numbers N j when the graphene layer is doped with a total filling factor ν. Assuming the broadening to be of Gaussian type with a width γ 0 the Fermi level is determined by solving the following equation:

ν = j Erf( E F -E j √ 2γ 0 ) (4)
from which one can calculate the individual filling factors

ν j = (1 + Erf( EF -Ej √ 2γ0
))/2 for each level E j . These E j will be used, later on, as fitting parameters dependent on the broken-symmetry phase under consideration. Note that in this approach all optical transitions to or from the n = 0 LL are allowed.

B. Description of the optical transitions

We first consider the transitions from the n = 0 LL to n = 1 LL. The Hamiltonian of the magneto-excitons, including their interaction with the Γ-phonon, denoted H (reminding that it describes the optical transitions allowed in the σ + polarization), is:

H =       ω 01 -E 1 0 0 0 g 1 √ ν 1 0 ω 01 -E 2 0 0 g 2 √ ν 2 0 0 ω 01 -E 3 0 g 3 √ ν 3 0 0 0 ω 01 -E 4 g 4 √ ν 4 g * 1 √ ν 1 g * 2 √ ν 2 g * 3 √ ν 3 g * 4 √ ν 4 ω ph       . ( 5
)
where ω 01 is the energy of the E 01 transition from the n = 0 LL to n = 1 LL in the absence of interactions and ω ph that of the Γ-phonon. This Hamiltonian describes the excitations from the 4 sublevels {E j , j = 1..4} of the n = 0 LL to the n = 1 LL. The matrix elements {g j , j = 1..4} respectively describe their interaction with the Γ-phonon, and is dependent on the wave-function character of the 4 sublevels (i.e., dependent on the broken-symmetry phase). In general, g j ∝ √ B [16], with a prefactor dependent on j and the brokensymmetry phase.

Similarly, the Hamiltonian describing the magnetoexcitons for the transitions from the n = -1 LL to the n = 0 LL (allowed in the σ -polarization) is written as:

H =       ω 01 + E 1 0 0 0 g 1 √ 1 -ν 1 0 ω 01 + E 2 0 0 g 2 √ 1 -ν 2 0 0 ω 01 + E 3 0 g 3 √ 1 -ν 3 0 0 0 ω 01 + E 4 g 4 √ 1 -ν 4 g * 1 √ 1 -ν 1 g * 2 √ 1 -ν 2 g * 3 √ 1 -ν 3 g * 4 √ 1 -ν 4 ω ph       . (6)
The total Hamiltonian H describing the magnetoexcitons is therefore:

H = H 0 0 H . (7)
We will also need to introduce the optical matrix elements M x and M y for corresponding transitions. These matrix elements depend on the ground state under consideration. In CDW case they are (see Sec. V A) :

M x /v 0 = { √ ν 1 , √ ν 2 , √ ν 3 , √ ν 4 , 0, - √ 1 -ν 1 , - √ 1 -ν 2 , - √ 1 -ν 3 , - √ 1 -ν 4 , 0}
and M y = iM x . For a different scenario, the optical matrix elements will be transformed to a different basis, as will be detailed in Sec. V. The Green's function for the magneto-excitons, is obtained as G = (( ω + iγ 01 ).I -H) -1 (where I is the unit matrix and γ 01 the broadening of the E 01 transition). This allows us to calculate the different components of the conductivity:

σ xx (ω) = i ω M T x .G.M x σ xy (ω) = 1 ω M T x .G.M * y (8)

V. OPTICAL CONDUCTIVITY IN THE DIFFERENT PHASES

Here, we calculate the optical conductivity for the different symmetry broken phases, using Eq. 8.

A. Charge density wave (CDW) phase

The CDW phase is characterized, at filling factor ν = 0, by two electronic LL full in one valley (say K ′ for instance) and two LL empty in the other valley (K). We will introduce a valley asymmetry ∆ V mainly determined by electron-electron interactions [9] and a Zeeman splitting ∆ S . Therefore the sequence of sublevels take the following form:

E 1 (K ′ ↑) = -∆ V /2 -∆ S /2 E 2 (K ′ ↓) = -∆ V /2 + ∆ S /2 E 3 (K ↑) = ∆ V /2 -∆ S /2 E 4 (K ↓) = ∆ V /2 + ∆ S /2 (9)
In this case, the parameters governing the electron-Γphonon interaction g 1 , g 2 on one hand and g 3 , g 4 on the other hand are of opposite sign. That is, where Ψ † K,s,n is the creation operator for electrons in valley K, spin s, Landau level n. On the other hand, the electron-light interaction (which determines M ) has the same sign at both valleys.

GS CDW + Γphonon|H e-ph Ψ † Ks,1 Ψ Ks,0 |GS CDW = g ph / √ 2 GS CDW + Γphonon|H e-ph Ψ † K ′ s,1 Ψ K ′ s,0 |GS CDW = -g ph / √ 2 (10
GS CDW + photon|H e-light Ψ † Ks,1 Ψ Ks,0 |GS CDW = 1 GS CDW + photon|H e-light Ψ † K ′ s,1 Ψ K ′ s,0 |GS CDW = 1 (11)
The results obtained for this phase are presented in Fig. 4, assuming ∆ V proportional to B, for two extreme values of the carrier density. The electron-phonon coupling was taken to be g ph = 2.3 × B[T ] meV, which agrees with density functional theory (DFT) calculations [23] and experiments [24][25][26]. The results are not very dependent on N s . The value of ∆ S = 0.15 meV B corresponds to a g-factor of 2.6 to be compared with 2.7 ± 0.2 reported in [22]. The splitting of the transition is directly governed by the amplitude of ∆ V whereas the introduction of ∆ S modifies only the relative amplitude of the two transitions. In all cases both ∆ V and g ph need to be finite to observe the effect. We finally note that, in this case, ∆ V > ∆ S in coherence with the assumption made in our previous work [14]. However there is no clear consensus about the field dependence on ∆ V [9]. Therefore one can alternatively assume that ∆ V is proportional to √ B. The corresponding results are displayed in Fig. 5 keeping all other parameters fixed. We obtained essentially the same results as in Fig. 4. Within the experimental errors we will not be able to differentiate between the two magnetic field variations of ∆ V .

The CDW state is compatible with the experimental results as we will see below.

B. Kekulé-distortion (KD) phase

In this phase [9], the K and K ′ valleys hybridize into linear combinations K, K ′ . At ν = 0, both spin ↑ and spin ↓ electrons occupy one of these valley-combinations, say K. The ν = 0 ground state for the KD phase is Ψ † K↑,0 Ψ † K↓,0 |0 . Therefore, the "natural" basis for this phase, where the density matrix is diagonal, is {K ↑, K ↓ , K ′ ↑, K ′ ↓} in contrast to the basis {K ↑, K ↓, K ′ ↑ , K ′ ↓} used in the CDW phase. Therefore the sequence of sublevels take the following form:

j = 1 : K ↑ j = 2 : K ↓ j = 3 : K ′ ↑ j = 4 : K ′ ↓ (12)
The transformation rules for the operators in this basis are

Ψ † K,s,n = 1 √ 2 Ψ † K,s,n + e iφ Ψ † K ′ ,s,n Ψ † K ′ ,s,n = 1 √ 2 Ψ † K,s,n -e iφ Ψ † K ′ ,s,n (13)
where Ψ † K,s,n is the creation operator for electrons in valley state K, spin s, Landau level n. Making use of this change of basis (Eq. 13) and Eq. 10 and Eq. 11, we derive that the electron-light matrix elements do not change with respect to the CDW phase, and the electron-phonon matrix elements g vanish by symmetry. That is, 14) and the same for K ′ . The structure of the Hamiltonian (Eq.7) becomes only diagonal and no splitting is observed when calculating the conductivity. Therefore the KD phase does not explain the experimental results.

GS KD + Γphonon|H e-ph Ψ † K↑,1 Ψ K↑,0 |GS KD = 0 GS KD + photon|H e-light Ψ † K↑,1 Ψ K↑,0 |GS KD = 1 (

C. Ferromagnetic (F) phase

In the F phase [9], the ground state, at filling factor ν = 0, is composed in both valleys K and K ′ of a single full LL with the same spin. In analogy with the CDW phase, we will introduce a valley asymmetry ∆ V and a Zeeman splitting ∆ S . Therefore the sequence of energy levels take the following form:

E 1 (K ′ ↓) = -∆ V /2 -∆ S /2 E 2 (K ↓) = ∆ V /2 -∆ S /2 E 3 (K ′ ↑) = -∆ V /2 + ∆ S /2 E 4 (K ↑) = ∆ V /2 + ∆ S /2 (15)
Note that, in this case, ∆ S should be larger than ∆ V to preserve the ferromagnetic nature of the state. In the present case the parameters governing the electron-Γphonon interaction (Eq.5,6) g 1 , g 3 on one hand and g 2 , g 4 in the other hand are of opposite sign.

The results are displayed in Fig. 6 where we have taken for ∆ S the same evolution that in the CDW phase and ∆ V ∝ B. The conductivity does not show any significant splitting of the main line. In fact there is an eigenvalue of the corresponding Hamiltonian larger than that of the main line but it remains optically inactive. Therefore here also, the F phase does not explain the experimental results.

D. Canted anti-ferromagnetic (CAF) phase

The CAF phase for the ground state is described by a spin in direction θ K in valley K and a spin in direction θ K ′ in valley K ′ . (The directions θ K and θ K ′ are in general not opposite to each other except in the special case of the anti-ferromagnetic phase). The direction θ K is oriented at an angle θ relative to the magnetic field B and the direction θ K ′ at an angle -θ with respect to it. (In the anti-ferromagnetic phase, θ = π/2). Here the Zeeman splitting should vary like ∆ S ∝ cos θ and if θ is close to π/2 this term should not play a dominant role. We choose the following order of states:

E 1 (K, θ K ) = -∆ 1 /2 -∆ 2 /2 E 2 (K ′ , θ K ′ ) = -∆ 1 /2 + ∆ 2 /2 E 3 (K, π + θ K ) = ∆ 1 /2 -∆ 2 /2 E 4 (K ′ , π + θ K ′ ) = ∆ 1 /2 + ∆ 2 /2 (16)
where the introduction of ∆ 1 reflects the CAF pattern of spin. We assume in addition that the asymmetry between valleys is reflected by ∆ 2 (favoring here the K valley). To preserve the CAF phase ∆ 2 should be smaller than ∆ 1 . Similar to the F phase, the parameters governing the electron-Γphonon interaction (Eq.5,6) g 1 , g 3 on one hand and g 2 , g 4 in the other hand are of opposite sign. The results are displayed in Fig. 7 where we have taken, ∆ 1 and ∆ 2 proportional to √ B. The results are not very dependent on the carrier concentration. We observe indeed a splitting of the transition when both ∆ 2 and g ph are different from zero : in fact the splitting is governed by ∆ 2 . In the present case we do not have, a priori, a guide for choosing the values of ∆ 1 and ∆ 2 . In order to be consistent with experimental results, we have taken for ∆ 1 a value which provides an upper transition energy close to that observed.

-3 I 1P2 (T) ∆ 1 (eV) Q 1.2 10 -3 R 1S2 (T) a) σ U V (a.W ) CAF phase X Y 35T Energy (eV) ` b 25T
However the evolution of the spectra does not reflect the experimental observations: whatever is the choice of parameters, the intensity of the high energy transition never reaches that of the main transition in contrast to the CDW phase where it should become dominant at fields higher than 35 T. This is discussed further in the next section.

VI. COMPARISON OF EXPERIMENT AND THEORY

The KD, F and CAF phases result in transmission spectra incompatible with experiment (Fig. 8). In the KD phase, electrons occupy linear combinations of the K and K ′ valleys; the electron-phonon matrix elements vanish by symmetry, resulting in a single transmission line. For the F phase, the occupancy of the K and K ′ valleys are almost equal (Sec. V C), and there is no significant splitting of the main transmission line (Fig. 8) . Similarly, the calculated transmission spectra for the CAF phase show a second CR line of much lower intensity than the main CR line. Deconvoluting these spectra with two Lorentzians we find a ratio of the CR weights of 0.9 ± 0.05 for the experiment, to be compared to the value 0.9 for the CDW phase and 0.4 for the CAF phase. Despite the introduction of valley asymmetry into the CAF phase, we find that it cannot explain the observed evolution of CR energies in the experiment.

The CDW phase has unequal occupation numbers of the n = 0 LL at the K and K ′ valleys, corresponding to a density modulation of the graphene A and B sublattices in real space. Unlike the ideal disorder-free CDW discussed in [9], both K and K ′ valleys have non-zero occupation number in our calculation, due to disorderinduced broadening. Nevertheless, the mechanism giving rise to the splitting of the E 01 transmission line remains essentially the same as illustrated by the simple model (Eq. 1) above.

The spin and valley splittings ∆ S and ∆ V parameterize our model for the CDW phase and determine the filling factors that enter the Hamiltonian and the optical matrix elements. We fix ∆ S using the experimental graphene gfactor measured in Ref. [22]. We treat ∆ V as a fitting parameter, obtaining ∆ V = 0.24 × B[T]meV. In our calculations, we use g ph = 2.3 × B[T]meV, which is in good agreement with density functional theory (DFT) calculations [23] and experiments [24][25][26]. We take the position of the E 01 transition line and its broadening to be given by v F = 1.01 × 10 6 ms -1 and γ 01 [meV] = 3 + 0.8 B[T] respectively, consistent with their measured values at low magnetic fields away from the Γ-phonon frequency. For the parameter γ 0 characterizing the broadening of the Landau levels in Eq.4, we have taken γ 0 = γ 01 /2 because the broadening of the E 01 transition should have contributions from both the n = 0 and n = 1 Landau levels. The calculations for the splitting at high mag-netic fields are in excellent agreement with the experimental transmission spectra for both samples S4 (Fig. 9) and S5 (Fig. 1). We neglect K-phonon absorption [14], which might account for the discrepancy between theory and experiment at the lower frequency and magnetic field range of our data. Upon changing the carrier density N s by a factor of ∼ 10 in our calculations, we find only minor changes in the transmission spectra, which is consistent with the observation that samples S4 and S5 have remarkably similar transmission spectra despite having different carrier densities. This is because the disorder induced broadening reduces the dependence of the sublevel filling factors (ν K↑ , etc.) on the carrier density.

For symmetry breaking driven by electron-electron interactions, the details of the screening function plays a vital role in determining the nature of the ground state. The additional screening afforded by the multiple graphene layers in our epitaxial graphene samples might favor the CDW configuration, with two electrons on the same sublattice, over the CAF state observed in hBNsupported samples [6]. Furthermore, coupling between rotationally misaligned layers breaks the local A-B sub-lattice (i.e. valley) symmetry [29,30], promoting the CDW ground state.

VII. CONCLUSIONS

conclusion, we have used magneto-optical spectroscopy to characterize a SU(4) symmetry broken phase in our epitaxial graphene samples. Based on the evolution of the transmission lines near the Γ phonon frequency, we identify this phase as a CDW phase for the specific samples considered, with different occupation numbers at valleys K and K ′ . Because of the valleysensitive nature of the electron-phonon interaction, the transmission study used here complements spin-sensitive transport measurements in tilted magnetic fields in the study of symmetry breaking in graphene. Our experimental method can be applied to open questions such as symmetry breaking of the different LLs in graphene and bilayer graphene, as well as the effect of disorder on the broken symmetry phase in these systems. where r,s scan the values 0 and ±1, 0 ≤ f r ≤ 1 is the occupation factor of the LL r, M r,s the optical matrix element, E r,s = E r -E s = E 01 and Γ rs (B) = γ 01 (B) measures the broadening of the transition. M r,s ∝ v 0 where v 0 is the Fermi velocity given by LDA calculations [31]. In the present work, we have taken for all samples v 0 = 0.85 × 10 6 ms -1 . This is a different parameter from v F which appears in E 01 because the energies and wavefunctions are corrected to different extents by the electron-electron interaction [31]. This approach requires the knowledge of the number of effective active layers as well as their carrier densities N s (ν = N s Φ 0 /B, Φ 0 being the flux quantum) which, in turn, implies some approximations.

The multi-layer dielectric model assumes that each graphene sheet is uniformly spread over the sample. This is a strong assumption, difficult to justify a priori and we have been lead to correct it by assuming a mean coverage which, in the present case for samples S4 and S5, has been determined to be about 70 per cent. We next evaluate the number N ef f for each sample. In the range of magnetic fields 12 to 17 T, the relative transmission spectra (Fig. 2, top panel) reaches values above 1 which depends on the number N ef f : we have therefore a guide to estimate this quantity. We estimate N ef f = 7 for samples S4 and S5.

The carrier density N s for each layer is determined in the following way: one knows that, for 2 < ν < 6, upon increasing B, the intensity of the E 01 absorption starts to increase, at the expense of the intensity of the E 12 transition (E 12 = E 2 -E 1 ). The intensity does not change with B for ν < 2. Therefore, the disappearance of the optical transition E 12 corresponds to ν = 2. Following the transmission spectra as a function of B, one can evaluate the carrier density N sm for each layer m. This is an iterative process which converges reasonably (within 20 per cent) but has to be done independently for each sample. The value of N s1 for the layer close to the SiC substrate can be set arbitrary to 5 to 6 10 12 cm -2 as given by transport data on samples grown under similar conditions: this layer indeed and the two following ones do not contribute to the transition E 01 in the present experiment. Finally, in the range of magnetic field larger than 27 T, where we focus our attention in this paper, the number of optically active layers (for optical transitions involving the n = 0 LL) ranges between 3 to 4 for samples S4 and S5 with carrier densities ranging from 0.5 to 12 ×10 11 cm -2 .

2 FIG. 1 :

 21 FIG. 1: (color online). Left panel: Evolution of the E01 transition for different values of the magnetic field beyond 27 T. Experimental transmission data for sample S5 (open circles) is compared with calculated transmission spectra (red lines), for different values of magnetic field, using the proposed CDW model. Right panel: 1 -(relative transmission) measured in the experiment and calculated (red lines) for the CDW phase, for magnetic fields B = 27, 29, 31, 33, 35 T. Deconvolution of the experimental spectra into two Lorentzians is shown in blue dashed and dotted lines.

FIG. 2 :

 2 FIG. 2: (color online). Top: Relative transmission spectra of sample S5, for different magnetic field values up to 35 T. Bottom: Evolution of the E01 transition for different values of the magnetic field between 27 T and 34 T, for samples S4 (open dots) and S5 (full lines).

  FIG. 3: (color online). Top: CR energies, resulting from the deconvolution of the experimental transmission traces, as a function of √ B for samples S4 (full dots) and S5 (open dots). Error bars are similar for both samples. The full lines are a linear fit of the data for the low energy transition at low (F2(B)) and high (F1(B)) field. Bottom: Comparison of the variation of experimental CR energies for sample S5 (open dots), with that calculated for the CDW phase (open squares), as a function of √ B. The size of the open squares mimics the relative oscillator strength of the optical transition.

N S 10 10 11 m - 2 FIG. 4 :

 24 FIG. 4: (color online). CDW phase with ∆V ∝ B: evolution of the σxx component of the conductivity for different values of the magnetic field between 25 T and 35 T for a) a carrier density Ns = 0.5 × 10 11 cm -2 and b) Ns = 1 × 10 12 cm -2 . In both cases g ph = 2.3 × B[T ] meV.

∆

  FIG. 5: (color online). CDW phase with ∆V ∝ √ B: evolution of the σxx component of the conductivity for different values of the magnetic field between 25 T and 35 T for a) a carrier density Ns = 0.5 × 10 11 cm -2 and b) Ns = 1 × 10 12 cm -2 . In both cases g ph = 2.3 × B[T ] meV.

FIG. 6 :

 6 FIG. 6: (color online). Ferromagnetic phase with ∆V ∝ B: evolution of the σxx component of the conductivity for different values of the magnetic field between 25 T and 35 T for a carrier density Ns = 5.0×10 11 cm -2 and g ph = 2.3 meV/T 1/2 .

- 2 FIG. 7 :

 27 FIG. 7: (color online). CAF phase with ∆1 ∝ √ B and ∆2 ∝ √ B: evolution of the σxx component of the conductivity for different values of the magnetic field between 25 T and 35 T for a) a carrier density Ns = 0.5 × 10 11 cm -2 and b) Ns = 1 × 10 12 cm -2 . In both cases g ph = 2.3 meV/T 1/2 .

FIG. 8 :

 8 FIG. 8: (color online). Real part of the conductivity measured in the experiment (open dots) and calculated with the CDW, CAF, F and KD phases (red continuous lines), for magnetic fields B = 31T (left panel) and 35T (right panel). The Deconvolution of the model into two Lorentzians is shown for the CDW and CAF phases (blue dashed and dotted lines).

FIG. 9 :

 9 FIG. 9: (color online). Fit of the experimental data (open dots) for sample S4 with the conductivity model derived for the CDW phase (continuous lines) using ∆V = 0.24 meV/T and g ph = 2.4 meV/T 1/2 .

Acknowledgements L.Z.T. and the theoretical analysis were supported by the Theory Program at the Lawrence Berkeley National Lab through the Director, Office of Science, Office of Basic Energy Sciences, Materials Sciences and Engineering Division, U.S. Department of Energy under Contract No. DE-AC02-05CH11231. Numerical simulations were supported in part by NSF Grant DMR10-1006184. Computational resources were provided by NSF through Tera-Grid resources at NICS and by DOE at Lawrence Berkeley National Laboratory's NERSC facility. We acknowledge the support of this work by the France-Berkeley Fund and the European Research Council (ERC-2012-AdG-320590-MOMB). * Corresponding author: sglouie@berkeley.edu.

Appendix

In general, the detailed analysis of magnetotransmission spectra requires the use of a multi-layer dielectric model including all layer dielectric properties of the sample. In particular, for each graphene sheet, one has to introduce the corresponding components of the optical conductivity tensor σ xx (ω) and σ xy (ω). Here, the x and y-axis lie in the plane of the sample. For instance σ xx (ω), in a one-electron approximation, for transitions involving the n = 0 LL, is written as:

M 2 r,s (f r (B) -f s (B)) ω -E r,s (B) + iΓ rs (B) (17)


Carlos Sánchez

Eduardo Bolea-Fernández

Marta Costas-Rodríguez

Charles- Philippe Lienemann

Jose-Luis Todolí

Frank Vanhaecke
email: frank.vanhaecke@ugent.be

Direct lead isotopic analysis of bioethanol by means of multi-collector ICP-mass spectrometry with a total consumption sample introduction system

Keywords: Lead isotope ratios, direct isotopic analysis, ICP-MS interface, total sample consumption system, bioethanol

A method has been developed for the direct (no sample pre-treatment and/or isolation of the target analyte from the sample matrix) lead isotopic analysis of bioethanol samples via multicollector inductively coupled plasma-mass spectrometry (MC-ICP-MS). A total consumption sample introduction system, the so-called high-temperature Torch-Integrated Sample Introduction System (hTISIS), equipped with a PFA micro-nebulizer and a heated small-volume spray chamber, has been used for (i) reducing the analyte concentration required for obtaining accurate and precise lead isotope ratio results and (ii) mitigating the effect of the ethanolwater ratio on the extent of mass bias.

The results obtained when using the hTISIS have been compared to those obtained with a more conventional sample introduction system, i.e. a micro-nebulizer mounted onto a cyclonic spray chamber at room temperature. The performance of both introduction systems has been assessed for two different plasma interfaces. The Pt standard sampling cone has been combined with an X-type or H-type skimmer cone, respectively. The sensitivity achieved with the hTISIS was between 3-and 7.5-fold higher, depending on the ethanol-water ratio, than that with the conventional sample introduction system, thus permitting accurate lead isotope ratios to be obtained at lower concentration levels without degradation of the precision. The external precisions, reported as 2 times the relative standard deviation (2RSD), for 207 Pb/ 206 Pb and 208 Pb/ 206 Pb were 0.007% and 0.008%, respectively, whereas the internal precision was 0.007% (2 RSD) for both isotope ratios.

The effects of ethanol content and the hTISIS temperature on the extent of mass bias have been evaluated for the four instrument setups (different sample introduction system/skimmer cone type combinations). The combination of (i) internal correction using NIST SRM 997thallium as an internal standard relying on Russell's law and (ii) external correction using NIST SRM 981 -lead, prepared in 75% ethanol, in a sample-standard bracketing (SSB) approach was used for mass bias correction. Although bioethanol samples may contain different amounts of water, the correction described above enabled adequate correction for mass bias in ethanolwater matrices with a water content of 0 to 40% and thus, also for actual bioethanol samples, when using the hTISIS operated at 125°C and an X-type skimmer cone.

The methodology has been successfully validated by means of lead isotopic analysis of bioethanol samples spiked with a lead standard previously characterized isotopically. Finally, as a proof of concept, actual bioethanol samples have been analyzed and significant differences in the lead isotope ratios have been observed.

Introduction

Biofuels are currently considered an interesting alternative energy source as they bring about a lower burden on human health and the environment than fossil fuels, while they could offer a solution for the increasing demand of energy in view of the limited petroleum stock. [1][START_REF] Köpke | Handbook of biofuels production: processes and technologies[END_REF][START_REF] Walker | Bioethanol: Science and technology of fuel alcohol[END_REF] Among the variety of biofuels available, bioethanol has emerged as one of the most promising products, as its use involves a reduction of the emission of greenhouse gases (GHG), such as CO 2 , CH 4 and N 2 O, up to 75% compared with fossil fuels. [1][START_REF] Köpke | Handbook of biofuels production: processes and technologies[END_REF][START_REF] Walker | Bioethanol: Science and technology of fuel alcohol[END_REF][START_REF] Sánchez | [END_REF][5] Bioethanol can be used in its pure form within modified spark-ignition engines or blended with petrol (ethanol-fuel). In the latter case, the bioethanol in the blend acts as an efficient octaneboosting agent, thereby substituting for toxic chemical additives, such as methyl tert-buthyl ether (MTBE). 1,[START_REF] Walker | Bioethanol: Science and technology of fuel alcohol[END_REF]6 As a result of the advantages of bioethanol compared to fossil fuels, its production and consumption have grown exponentially during the last two decades. [START_REF] Lemos | Global Bioethanol[END_REF] Different biofuel generations have appeared along this time. The first-generation bioethanol is the alcoholic product generated from sugars, seeds or starch by microorganisms. 6 However, the fuel-food competition is regarded an important ethical issue impeding the global implementation of this first-generation bioethanol. 1,6,[START_REF] Fargione | [END_REF] The second-generation bioethanol overcomes this problem, since the raw material used is based on whole plant biomass or on non-edible crop production. [START_REF] Köpke | Handbook of biofuels production: processes and technologies[END_REF]6,9 However, the production process is complex and the yield reached is lower than that of the first-generation processes. 6 The final product generated from these processes contains mainly ethanol and a variable percentage of water. It may also contain organic and inorganic pollutants (including metals) at minor and trace levels. Nowadays, only limited regulations for the levels of metals in bioethanol samples exist [START_REF]ASTM International, ASTM D4806-16a, Standard Specification for Denatured Fuel Ethanol for Blending with Gasolines for Use as Automotive Spark-Ignition Engine Fuel[END_REF][START_REF]ABNT NBR 11331-Ethyl alcohol -Determination of iron and copper concentrations -Atomic absorption spectrophotometry method[END_REF] , although recent works paid attention to the analysis of this type of samples via inductively coupled plasma-optical emission spectrometry (ICP-OES) [START_REF] Sánchez | [END_REF][START_REF] Sánchez | [END_REF] or inductively coupled plasma-mass spectrometry (ICP-MS) [START_REF] Sánchez | [END_REF]13 . However, not only elemental analysis of bioethanol is relevant, e.g., for the purpose of quality control, but also isotopic analysis can provide valuable information, e.g., related to the geographical provenance and type of raw material used for its production.

In general, the conventional methodology used for high-precision isotopic analysis via multicollector inductively coupled plasma-mass spectrometry (MC-ICP-MS) involves the isolation of the analyte from the matrix using off-line chromatographic methods. [14][15][16] This approach reduces the occurrence of spectral interference and allows a reliable correction for the mass bias caused by instrumental mass discrimination. However, this procedure is laborious, timeconsuming and costly, and it can generate a considerable amount of waste. As an alternative, some papers describe the use of high-performance liquid chromatography (HPLC) coupled to MC-ICP-MS for on-line chromatographic separation. 17,18 Even after analyte isolation, adequate correction for mass bias often involves a combination of internal correction relying on the use of an admixed internal standard and external correction in a sample-standard bracketing approach (SSB). 16,[19][20][21][22] High-precision isotopic analysis of Pb by multi-collector ICP-MS is used for provenancing (determination of the geographical origin) purposes. [22][23][START_REF] Reimann | Hoogewerff and The GEMAS Project Team[END_REF][25] Lead isotope ratios have been determined in fossil fuels and other organic samples 26,27 using the conventional methodology, involving off-line chromatographic isolation of Pb from the sample matrix. However, to the best of the authors' knowledge, so far, isotopic analysis has not been performed for biofuels in general, and bioethanol in particular.

This study explored Pb isotopic analysis of bioethanol, whereby an effort was done to avoid extensive sample pretreatment. The direct introduction of organic samples in a MC-ICP-MS instrument is challenging due to several reasons: the organic matrix can negatively affect the stability of the plasma, while the organic solvents can exert an effect on the extent of mass bias that might not be adequately corrected for due to the difference in matrix composition between the external standard and the sample, thus possibly jeopardizing the accuracy of the results obtained.

The total consumption sample introduction system hTISIS (high-temperature Torch-Integrated Sample Introduction System) 28 was shown to be able to overcome the main problems related with elemental analysis of organic samples in general, and of fuels and biofuels in particular. The hTISIS operated under optimized conditions can minimize matrix effects caused by ethanol or other organic solvents delivered to the plasma, because the differences in transport efficiency are mitigated or even eliminated. [START_REF] Sánchez | [END_REF]13,29,30 In addition, the use of the hTISIS contributes to an enhanced plasma stability, since the optimum sample uptake rates used with the hTISIS are lower than those used with the conventional sample introduction systems and, thus, also the volume of sample required for analysis is reduced. Moreover, the sensitivity provided by the hTISIS is generally higher than that with the conventional sample introduction systems. The combination of an hTISIS and MC-ICP-MS has been previously used for carrying out studies in aqueous sample matrices, aiming at improving the sensitivity and/or at reducing the volume of sample required for the analysis. [31][32][33][34][35] However, its application to direct isotopic analysis of organic samples has not been investigated yet.

The main objective of this work was to explore the analytical capabilities of the hTISIS as a sample introduction system in the direct, i.e. without sample pre-treatment and/or chromatographic separation, isotopic analysis of lead in bioethanol samples via MC-ICP-MS. For this purpose, the effect of various ethanol-water mixtures on the lead isotope ratio results was evaluated for the hTISIS (operated at different temperatures) and a conventional sample introduction system (microconcentric nebulizer and cyclonic spray chamber) coupled to a MC-ICP-MS unit using two different skimmer cone types (X and H -type).

Experimental 2.1. Aqueous standards and certified reference materials

Isotopic reference materials of lead (NIST SRM 981) and thallium (NIST SRM 997) from the National Institute for Standards and Technology (NIST, Gaithersburg, MD, USA), were used for the correction of the mass bias induced by instrumental mass discrimination; they will be referred to along the manuscript as "NIST-Pb" and "NIST-Tl", respectively. The isotope ratios reported in the certificate, and thus used for mass bias correction are 208 Pb/ 206 Pb = 2.1681 ± 0.0008, 207 Pb/ 206 Pb = 0.91464 ± 0.00033, 208 Pb/ 207 Pb = 2.3704 ± 0.0009 and 205 Tl/ 203 Tl = 2.38714 ± 0.00101.

An in-house standard solution of lead (Inorganic Ventures, Christiansburg, VA, USA, lot D2-PB03020) was characterized for its isotopic composition using a conventional approach, i.e. by means of MC-ICP-MS with a conventional sample introduction system relied on Tl for the internal standardization using NIST-Tl and SSB approach using NIST-Pb as standards 36 . This standard solution will be termed "IH-Pb". The lead isotope ratios for IH-Pb are 208 Pb/ 206 Pb = 2.17616 ± 0.00019, 207 Pb/ 206 Pb = 0.92331 ± 0.00006 (average ± 95% confidence level for n=15 and α=0.05, determined in a 3% HNO 3 matrix).

Ethanol-water standards and bioethanol samples

Different ethanol-water standards containing ethanol concentrations ranging from 10 to 100% (v/v) were prepared using ultrapure water (resistivity ≥ 18.2 MΩ cm) obtained from a Direct-Q water purification system (Millipore, Molsheim, France) and absolute ethanol of analytical reagent grade (Fisher Scientific, Waltham, MA, USA). Isotopic reference materials NIST-Pb and NIST-Tl were added into the ethanol-water solutions in a 1:1 ratio. Additionally, standards containing IH-Pb and NIST-Tl were prepared in 50%, 75% and 100% of ethanol. The Pb and Tl concentrations added to each ethanolic solution were adapted such as to obtain the same signal intensity in every ethanol-water mixture.

Sample-standard bracketing solutions containing NIST-Pb and NIST-Tl in a 1:1 ratio were prepared in 75% ethanol. The Pb:Tl concentration ratio was selected aiming at obtaining signal intensities in the same order of magnitude for Pb and Tl isotopes, as it is required for an adequate mass bias correction. The Pb concentrations at which the isotope ratio measurements were carried out were: (i) 10 µg L -1 with the hTISIS and X-type skimmer cone; (ii) 15 µg L -1 for the conventional introduction system and the X-type skimmer cone; (iii) 40 µg L -1 for the hTISIS and the H-type skimmer cone; and (iv) 60 µg L -1 for the conventional introduction system and the H-type skimmer cone. The concentrations were different for different sample introduction systems and skimmer cones in order to reach the same signal intensity for the four approaches.

For method validation purposes, seven bioethanol samples containing lead concentrations lower than 0.025 µg L -1 were spiked with 5 µg L -1 of IH-Pb and 5 µg L -1 of NIST-Tl. Furthermore, lead isotope ratios were measured in six real bioethanol samples. The raw material used for the manufacturing of bioethanol can be of different nature; bioethanol is mainly obtained from sugarcane, wheat, wine alcohol, or biomass (when produced from the latter, it is termed second-generation bioethanol). Some of these samples were dehydrated bioethanol, whereas some of them contained water in a content ranging from approximately 2 to 35%. Tl concentration in the samples was below 0.003 µg L -1 (method detection limit for thallium quantification-MDL), see section 2.3.

Instrumentation and measurements

An Agilent 7900 ICP-MS Spectrometer (Agilent, Santa Clara, CA, USA) was used for Pb and Tl quantification in the bioethanol samples. Lead and thallium standards for external calibration were prepared in 2% HNO 3 by appropriate dilution of a stock solution containing 1 g L -1 of Pb (Inorganic ventures) and 1 g L -1 of Tl (Inorganic ventures). Intensities for 203 Tl, 205 Tl, 206 Pb, 207 Pb and 208 Pb were monitored in vented mode (without addition of collision and/or reaction gas).

Pb and Tl isotope ratios were measured using a Thermo Scientific (Bremen, Germany) Neptune MC-ICP-MS instrument equipped with a 130 m 3 /h dry interface pump, operated at low mass resolution. All solutions were delivered to the nebulizer by means of a peristaltic pump (Spetec, Erding, Germany) and by using a 0.25-mm id flared-end solvent-resistant tubing (Glass Expansion, Melbourne, Australia).

Samples and standards were introduced into the plasma by means of a PFA nebulizer (Elemental Scientific, Omaha, NE, USA) coupled onto a spray chamber. Two different sample introduction systems were evaluated: (i) an hTISIS with a 9 cm 3 single-pass spray chamber and operated at different temperatures [START_REF] Sánchez | [END_REF]13,28,29,37 , and (ii) a cyclonic spray chamber (Elemental Scientific, Omaha, NE, USA) operated at room temperature, as a conventional sample introduction system. The instrument settings and the data acquisition parameters selected are gathered in Table 1. For final method validation and for analysis of the real bioethanol samples, the hTISIS was operated at 125°C. Platinum sampling and skimmer cones were used because of the introduction of O 2 for adequate handling of the ethanolic solutions. Two different skimmer cone types were tested, the H-type (conventional model) and X-type (providing higher sensitivity). 38 A quartz torch with a 1.0 mm inner diameter injector was used. A Smart Mass Flow 5850S (Brooks Instruments, Hatfield, PA, USA) mass flow controller, connected to a Brooks Microprocessor Control & Read Out Unit 0152 (Brooks Instruments, Hatfield, PA, USA), was used for introducing 0.2 L min -1 of additional gas (20% of O 2 + 80% of Ar) to avoid carbon deposition on the torch and the interface cones. [START_REF] Sánchez | [END_REF]13,39 Correction for the bias caused by instrumental mass discrimination was accomplished using the combination of internal correction relying on Tl as an internal standard, followed by external correction in a sample-standard bracketing (SSB) approach. For the internal correction with Tl as an internal standard, different correction approaches were evaluated: (i) Russell's law 16 , (ii) the method proposed by Woodhead et al. 19 , and (iii) the latter method as further refined by Baxter et al. 20 An ethanol-water mixture with 75% ethanol was consistently used as matrix for the SSB standard. In all cases, the Pb signals for the blanks were negligible compared to those for samples and standards (<0.1%).

Results and discussion

Effect of sample introduction system and skimmer type on the sensitivity

The use of different sample introduction systems and/or plasma-spectrometer interfaces strongly affects the sensitivity of ICP-MS instruments. In this section, the effects of (i) the Pt skimmer cone type (X vs H) and (ii) the sample introduction system (conventional vs hTISIS) on the absolute sensitivity of 208 Pb in 10, 50 and 100% ethanol have been evaluated. As can be seen in Figure 1.a., a 5 to 7-fold improvement in sensitivity was obtained when replacing the H-type by the X-type skimmer cone. This enhancement in sensitivity is caused by the different geometry of the X skimmer cone, which provides a significantly higher ion extraction efficiency from the plasma. 40,41 Figure 1.b. shows the sensitivity obtained for different ethanol-water matrices using the hTISIS operated at 125°C (the optimized temperature, providing maximum sensitivity) and the conventional introduction system using the X-type skimmer cone. It is noteworthy that the use of the hTISIS leads to an enhanced Pb signal sensitivity for all of the ethanol-water standards tested. A 1.6-fold increase was obtained at 100% ethanol, whereas an improvement factor of 2.3 and 3.7 was found for the Pb standards in 50 and 10% of ethanol, respectively. These differences in the enhancement factor can be attributed to the different operation mechanism of the sample introduction systems used. In the case of the hTISIS working under optimum conditions, the analyte transport efficiency is virtually 100% for any ethanol content. However, the transport efficiency with the conventional sample introduction systems strongly depends on the physical properties of the sample as these define the nebulization efficiency and aerosol droplet size distribution, with the sensitivity increasing with increasing ethanol concentration. It needs to be stressed that the hTISIS provided the aforementioned improvement despite working at half of the sample uptake rate than that of the conventional introduction system. Taking this fact into account, the improvement in sensitivity was 3.2-, 4.6and 7.4-fold for Pb standards in 100, 50 and 10% of ethanol, respectively. The same trend in the sensitivity was observed for Tl.

Next to this enhancement in sensitivity, with the hTISIS the sensitivity also depends much less on the ethanol concentration of the samples, which is advantageous for the isotopic analysis of bioethanol samples, since these may be characterized by different water contents. [START_REF] Sánchez | [END_REF]13 Under optimum conditions, the sensitivity was 0.07-0.08 V per µg L -1 (0.14-0.16 V ng -1 ) for 208 Pb and 0.10-0.11 V per µg L -1 (0.20-0.22 V ng -1 ) for 205 Tl with amplifiers equipped with a 10 12 Ω resistor. As is already known, a reduction of analyte concentration degrades the isotope ratio precision.

Figure 2 shows the influence of the Pb concentration on the 208 Pb/ 206 Pb ratio in standard solutions prepared in 100% ethanol with the hTISIS (at 125°C) and with the conventional sample introduction system, both with the X-type skimmer cone (higher sensitivity without compromising precision compared to H-type skimmer cone). The accuracy and precision of the isotope ratio at different lead concentrations can be estimated from the average value and the relative standard deviation for three measurement replicates of NIST-Pb in an SSB approach. In terms of precision, the hTISIS provided better results than did the conventional sample introduction system at low Pb concentration (see Figure 2), due to the sensitivity enhancement provided by this total sample consumption system. Additionally, the standard deviation was similar at all concentrations > 2 µg L -1 with the hTISIS, whereas it increased significantly at concentrations < 10 µg L -1 in case of the conventional sample introduction system. Based on these results, it can be concluded that the lowest Pb concentration required for obtaining precise lead isotope ratios in 100% ethanol matrix was 2 µg L -1 for the hTISIS operated at 125°C and 10 µg L -1 for the conventional sample introduction system.

When using the H-type skimmer cone, the minimum concentration of Pb required to obtain precise lead isotope ratios in a 100% ethanol matrix is about 15 µg L -1 for the hTISIS and and 25 µg L -1 for the conventional sample introduction system. Moreover, for both types of skimmer cone, the minimum Pb concentration required remains constant for any ethanol-water ratio in case of hTISIS, whereas it increases with the water content in case of the conventional sample introduction system (see Figure 1).

Effect of sample introduction system and skimmer type on the mass bias correction.

Different approaches for correction for instrumental mass discrimination were evaluated for different ethanol-water matrices, always using a standard in 75% ethanol for external correction of any ethanol-water mixture, as detailed in section 2.3. The different mass-bias correction methods, all based on internal correction relying on the use of Tl as an internal standard followed by external correction in an SSB approach did not provide results showing statistical differences and, finally, the Russell's law (followed by SSB) was relied on for correction of instrumental mass discrimination. 16 Figure 3 provides corrected lead isotope ratios obtained for all of the ethanolic matrices (concentrations between 10 and 100% of ethanol), and for the four combinations of sample introduction system/interface (Figure 3).

The mass spectrometer interface plays an important role in the phenomenon of mass discrimination. 41 With the H-type skimmer cone, mass bias correction, relying on Russell's law (NIST-Tl) followed by SSB (NIST-Pb prepared in 75% ethanol), was found to be adequate for ethanol-water mixtures containing from 10 to 100% and from 30 to 100% of ethanol for the conventional sample introduction system and the hTISIS, respectively. However, as indicated in previous sections, the sensitivity is strongly deteriorated when using the H-type skimmer cone, which means that this setup is less suited for lead isotopic analysis of bioethanol samples, since the concentration of the target analyte is generally lower than 3 µg L -1 . [START_REF] Sánchez | [END_REF] As shown in section 3.2, the minimum concentration required for obtaining accurate and precise lead isotope ratios using this skimmer cone type is 15 µg L -1 for the hTISIS and 25 µg L -1 for the conventional sample introduction system, respectively. Thus, the use of the X-type skimmer cone is advisable for obtaining reliable lead isotope ratio measurements at low concentration levels (2 µg L -1 of Pb). However, it was found that the combination of X-type skimmer cone and cyclonic spray chamber did not enable adequate correction for mass bias, unless in the case of a close matrix-matching, i.e., when the ethanol content of standards and samples is approximately the same (see Figure 3). The Pb isotope ratios obtained for almost all of the ethanol-water mixtures were statistically different from the certified value for NIST-Pb. Only for samples with an ethanol content between 70 and 90%, the 208 Pb/ 206 Pb ratio was not significantly different from the certified value. It has been suggested before that the X-type skimmer cone geometry is associated with additional contributions to the instrumental mass fractionation, thus requiring a closer concentrationmatching of sample and standard solutions. 41 It could therefore be concluded that the influence the ethanol content exerts on the mass bias prevents adequate correction for mass bias using internal correction relying on admixed Tl standard (NIST-Tl), followed by the SSB approach with a NIST-Pb in 75% ethanol when using the conventional sample introduction system.

In contrast, when using the hTISIS at 125°C combined with the X-type skimmer, the effect of the ethanol content on the extent of mass discrimination did not prevent adequate correction in a range from 60 to 100% of ethanol. It should be noted that for bioethanol samples after the first distillation step, a water content higher than 40% is not expected. Thus, based on (i) the improvement in sensitivity and (ii) the possibility to adequately correct for the mass bias in the range of 60 to 100% of ethanol, the combination of hTISIS and X-type skimmer was selected as the most suitable setup to accomplish Pb isotopic analysis of bioethanol samples.

Effect of hTISIS temperature on the mass bias correction

The temperature at which the hTISIS is operated is a relevant parameter, as has been previously observed for elemental [START_REF] Sánchez | [END_REF]13,29,30 and isotopic analyses. [31][32][33] The hTISIS temperature was optimized for ethanol-water mixtures representative for actual bioethanol samples, i.e., with an ethanol content from 60 to 100% while using the X-type skimmer. The 208 Pb/ 206 Pb ratios obtained upon combined internal (NIST-Tl) and external (NIST-Pb in 75% ethanol) correction are shown in Figure 4 It can be observed that, within the ethanol concentration range evaluated (60 -100%), the mass bias is not completely corrected for, neither with the conventional sample introduction system (Figure 4.a), nor with the hTISIS at temperatures lower than 125°C (Figure 4.b and Figure 4.c). This can be explained by differences in transport efficiency for different ethanolwater matrices. This hypothesis is in agreement with the partial mitigation of the effect of the ethanol concentration on the extent of mass bias as observed when the hTISIS is heated from room temperature (Figure 4.b) to 75°C (Figure 4.c). In addition, the Pb isotope ratio precision is deteriorated with the hTISIS in comparison with the conventional sample introduction system, when the hTISIS is operated at temperatures lower than 125°C, which can be related to sudden evaporation of ethanol droplets impacting on the chamber walls at lower temperatures.

As has been demonstrated (see section 3.3), when the hTISIS was heated to 125°C, the mass bias could be adequately corrected for in the range of interest (60 -100% of ethanol) using a single Pb standard in 75% ethanol for external correction and good precision was achieved (Figure 4.d). Under these conditions, the transport efficiency is virtually the same, independent of the ethanol concentration.

Finally, an increase of the hTISIS temperature above 125°C (Figure 4.e and Figure 4.f) resulted in a degradation of lead isotope ratio accuracy and, especially, in isotope ratio precision. It is hypothesized that further heating of the aerosol does not affect the transport efficiency, but is having a negative effect on the plasma stability thus giving rise to an increase in the signal noise.

All these observations are in agreement with the trends in sensitivity. The intensity was matrixdependent when the hTISIS was operated at temperatures lower than 125°C, being this effect more pronounced at room temperature than at 75°C. However, the sensitivity obtained for hTISIS temperatures ≥ 125°C was found to be the same for standards containing between 60% and 100% of ethanol. In addition, the optimum sensitivity was obtained at 125°C, reaching a plateau at this temperature.

Robustness of the method to real matrices

For method validation, seven bioethanol samples produced from four different raw materials and containing different water contents ranging from 0 to 30% (~100 to ~70% ethanol) and a commercial ethanol sample (96%, analytical grade), all containing < 0.025 µg L -1 of Pb (method detection limit for lead quantification -MDL) were used. It needs to be further stressed that the water content in bioethanol samples is unknown under normal conditions. The samples were spiked with 5 µg L -1 of IH-Pb standard solution and subsequently measured for their Pb isotopic composition using MC-ICP-MS with the hTISIS operated at 125°C and the X-type skimmer cone. As Pb was present in the bioethanol samples selected at ultra-trace concentration levels (<MDL), no "blank" subtraction was performed because its effect on the lead isotope ratios was considered negligible within the precision attainable.

The results were plotted versus the reference value for the IH-Pb standard solution, previously characterized isotopically ( 208 Pb/ 206 Pb = 2.17616 ± 0.00019, 207 Pb/ 206 Pb = 0.92331 ± 0.00006).

Figure 5 shows the 208 Pb/ 206 Pb ratio for the eight matrices spiked with 5 µg L -1 of the IH-Pb standard solution. The Pb isotope ratio data obtained were found to be in good agreement with the reference value and no significant differences were observed, neither among the individual results, nor between any of these results and the certified value. Thus, the use of the hTISIS operated at 125°C for sample introduction into the MC-ICP-MS unit equipped with the X-type skimmer cone, is a reliable approach for the direct isotopic analysis of lead in bioethanol samples.

Lead isotope ratios in bioethanol

As a proof of concept, 208 Pb/ 206 Pb and 207 Pb/ 206 Pb isotope ratios were determined in six real bioethanol samples manufactured from diverse raw materials and containing different water contents (from approximately 0 to 30% of water). The isotope ratio measurements were performed directly in the bioethanol samples without any sample pre-treatment i.e., without sample digestion and/or analyte isolation, and by using the hTISIS operated under optimum conditions (125°C) and with the X-type skimmer cone.

Pb concentrations in the bioethanol samples ranged from 0.8 to 2.5 µg L -1 . For three of the samples, the Pb concentration was lower than 2 µg L -1 , and thus, these samples were preconcentrated by evaporation at room temperature for reaching the minimum concentration required to carry out precise and accurate lead isotope ratios measurement (see section 3.2).

It must be ensured that the ethanol content is higher than 60% after the preconcentration step.

In Figure 6, the data obtained for the bioethanol samples are presented in a three-isotopes plot ( 207 Pb/ 206 Pb vs 208 Pb/ 206 Pb). All of the bioethanol samples showed a lighter lead isotopic composition than did the NIST-Pb and IH-Pb standards. 208 Pb/ 206 Pb and 207 Pb/ 206 Pb ratios in the samples ranged from 2.1024 ± 0.0013 to 2.1389 ± 0.0008 and from 0.8635 ± 0.0010 to 0.8927 ± 0.0006, respectively. These ratios are similar to those reported in the literature for wine and other alcoholic beverages coming from fermentation and distillation processes. 42,43 The second-generation bioethanol sample showed a heavier Pb isotopic composition than the other samples (first generation biofuel). Additionally, differences were observed between samples coming from different raw materials.

Several reasons could explain the results obtained: (i) Pb isotope ratios are governed by the corresponding ratios in the raw material and consequently, by the lead isotopic composition of the soil where this material grows and/or by that of airborne deposition [START_REF] Reimann | Hoogewerff and The GEMAS Project Team[END_REF][44][45][46] . (ii) It should also be verified whether the Pb isotopic composition is modified during some steps of the biomass processing to convert it in bioethanol (i.e., syrups extraction and purification, fermentation and distillation), as already demonstrated for metallurgical processes involving heating steps at moderate temperatures 47,48 . (iii) Finally, it should also be taken into account that the lead in the final bioethanol product could be a combination of different lead sources.

Recent studies suggested however that the main source of lead in bioethanol is the raw material and that the samples are not polluted during the production process. 49 In any case, further work is required to elucidate the sources of the lead isotopic variability in bioethanol.

However, these first results suggest that Pb isotopic analysis could be valuable in the field of biofuels for traceability purposes and for establishing the raw material used for their preparation and/or their provenance. In this work, the geographic origin of the samples was unknown and thus, further research is required for evaluating the potential of Pb isotopic analysis for provenancing purposes in this context.

Conclusions

A novel analytical method has been developed for the direct Pb isotopic analysis of bioethanol samples via hTISIS-MC-ICP-MS, without laborious sample preparation and/or analyte isolation.

The hTISIS showed to be advantageous compared to the conventional sample introduction system. The use of hTISIS operated at 125°C in combination with the use of a Pt X-type skimmer cone was found to offer the best conditions for Pb isotopic analysis of bioethanol samples with different ethanol concentrations.

Additionally, a close matrix matching was demonstrated not to be necessary as the instrumental mass bias can be adequately corrected for by using a combination of internal correction using a Tl internal standard, relying on the Russell's law, followed by external correction using a Pb standard in 75% ethanol in an SSB approach.

This approach is suitable for direct Pb isotopic analysis of bioethanol samples with a water content up to 40%. It provides accurate Pb isotope ratios at concentrations > 2 µg L -1 with a precision that is similar to those reported in previous papers (use of conventional sample introduction system). The isotope ratios obtained for the spiked bioethanol samples agreed well with the reference value.

As a proof of concept, Pb isotopic analysis of bioethanol samples coming from different raw materials and of unknown geographical origin revealed significant differences. Thus, the combination of the hTISIS with the MC-ICP-MS opens up new applications for the direct isotopic analysis of biofuel samples for traceability purposes and for obtaining additional

Figure 1 .

 1 Figure 1. Effect of ICP-MS interface (a) and introduction system (b) on the sensitivity. Volume hTISIS: 500 µL, volume conventional spray chamber: 1 mL. Error bars = 2s.

Figure 2 .

 2 Figure 2. Effect of lead concentration on accuracy and precision (H-type skimmer). (a) 208 Pb/ 206 Pb ratios for NIST-Pb. Error bars represent ± 2s (n=3). (b) RSD (%) for the 208 Pb/ 206 Pb ratio results plotted in Figure 2.a. Matrix: 100% ethanol.

Figure 3 .

 3 Figure 3. Effect of the matrix composition on the effectiveness of mass bias correction via a combination of internal correction (based on admixed Tl) and external correction using a Pb standard solution in 75% ethanol for the different skimmer types and sample introduction systems under optimum conditions. Error bars represent ± 2s (n=3). The shaded area shows the certified ratio ± confidence level.

  for the different temperatures studied: room temperature (Figure 4.b), 75°C (Figure 4.c), 125°C (Figure 4.d), 200°C (Figure 4.e) and 250°C (Figure 4.f). The results were also compared with those obtained with the conventional sample introduction system working at room temperature (Figure 4.a).

Figure 4 .

 4 Figure 4. Effect of hTISIS temperature on the effectiveness of mass bias correction. (a) Cyclonic spray chamber at RT (Ref. system), (b)-(f) h-TISIS at (b) RT, (c) 75°C, (d) 125°C, (e) 200°C, and (f) 250°C. Error bars represent ± 2s (n=3). The shaded area represents the certified ratio ± confidence level.

Figure 5 .

 5 Figure 5. 208 Pb/ 206 Pb ratio obtained for spiked bioethanol and ethanol samples with 5 µg L -1 of IH-Pb. Error bars represent ± 2s (n=3). The shaded area indicates the reference value ± confidence level.

Figure 6 .

 6 Figure 6. Three-isotopes plot for bioethanol samples coming from different raw materials. Error bars represent ± 2s.

Table 1 .

 1 Conditions used for isotope ratio measurements (Thermo Scientific Neptune MC-ICP-MS)

	Instrument settings
	Introduction system	hTISIS †	Conventional
	Sample uptake rate (µL min -1 )	50	100
	Spray chamber temperature (°C)	125	RT
	RF power (W)	1400	1300
	Plasma gas flow rate (L min -1 )	15
	Auxiliary gas flow rate (L min -1 )	0.80

† Optimum conditions in terms of sensitivity selected for the analysis of bioethanol samples.

RT: room temperature

Acknowledgements

The authors wish to acknowledge IFPEN and UNGDA for the samples provided. Carlos Sánchez thanks the Ministry of Education, Culture and Sports, Spain for the contract FPU 13/01438 and the mobility grant EST16/00433 and IFPEN for financial support. Marta Costas-Rodríguez thanks FWO-Vlaanderen for her postdoctoral grant. Eduardo Bolea-Fernandez thanks BOF-UGent for his postdoctoral grant.

Effect of sample introduction system and skimmer type on the isotope ratio precision and accuracy

The isotope ratio precision is a crucial parameter for isotopic analysis, since it strongly affects the overall uncertainty associated to each measurement result, and thus, the capability to reveal small variations in the isotopic composition, required for meaningful interpretation of the results.

Before assessing the effect of sample introduction system and skimmer cone type on the isotope ratio precision and accuracy for a given ethanol-water mixture (75% ethanol), further used as the SSB standard, the possible influence of the ethanol content on the external precision was evaluated by measuring the 208 Pb/ 206 Pb and 207 Pb/ 206 Pb ratios in NIST-Pb standards prepared in matrices containing ethanol concentrations in the range of 10 to 100%. For such experiment, the concentration of Pb for the hTISIS was set at 50 µg L -1 , while for the conventional sample introduction system, the Pb content was adapted such that a similar signal intensity was observed. No differences (or trend) of isotope ratios precision were observed as a function of the ethanol concentration, with isotope ratio precisions (reported as two times the relative standard deviation -2RSD) for each standard solution distributed randomly from 0.003 to 0.007% and from 0.001 to 0.006% for 208 Pb/ 206 Pb and 207 Pb/ 206 Pb, respectively (n=5, 2RSD).

Internal and external precisions (Table 2) were evaluated for the four instrumental setups: (i) hTISIS and X-type skimmer cone; (ii) hTISIS and H-type skimmer cone, (iii) cyclonic spray chamber and X-type skimmer cone and, (iv) cyclonic spray chamber and H-type skimmer cone. Internal precision was evaluated by using the relative standard deviation (2RSD) obtained for 42 consecutive measurements of a Pb standard solution prepared in 75% ethanol. The concentration of NIST-Pb was adapted for obtaining the same signal intensity with the four instrumental setups used (see Table 2). Significant differences in the precision of 208 Pb/ 206 Pb and 207 Pb/ 206 Pb ratios were not observed and similar precisions (2RSDs) were obtained for the four approaches, with values ranging from 0.007% to 0.010% for 208 Pb/ 206 Pb and around 0.007% for 207 Pb/ 206 Pb.

External precision was calculated as the relative standard deviation (2RSD) obtained for all measurements of the NIST-Pb standard (prepared in 75% ethanol) carried out during a complete measurement session (8h -12h) for each of the four instrumental approaches. The findings for external precision were similar to those observed for internal precision and no significant differences were noticed between the sample introduction systems and/or skimmer cone types used. Moreover, precisions obtained for 205 Tl/ 203 Tl and lead isotope ratios were not statistically different.