Fabrice Gerbier
email: fabrice.gerbier@lkb.ens.fr

Jean Dalibard
email: jean.dalibard@lkb.ens.fr

Gauge fields for ultracold atoms in optical superlattices

)], laser-assisted tunneling between adjacent sites creates an effective magnetic field. In the tight-binding approximation, the atomic motion is described by the Harper Hamiltonian, with a flux across each lattice plaquette that can realistically take any value between 0 and π. We show how to take advantage of the superlattice to ensure that each plaquette acquires the same phase, thus simulating a uniform magnetic field. We discuss the observable consequences of the artificial gauge field on non-interacting bosonic and fermionic gases. We also outline how the scheme can be generalized to non-Abelian gauge fields.

Introduction

The fractional quantum Hall (FQH) phases realized by two-dimensional electron gases in very large magnetic fields are among the most intriguing states of matter (see, for instance, [START_REF] Girvin | The quantum Hall effect: Novel excitations and broken symmetries[END_REF]). In such systems, electrons "bind" to magnetic vortices, forming strongly correlated phases with striking properties, such as exotic excitations ("anyons") which obey fractional statistics [START_REF] Stern | Anyons and the quantum Hall effect -a pedagogical review[END_REF]. Analogous quantum Hall phases should also arise in cold atomic gases when they are set into fast rotation (see [START_REF] Cooper | Rapidly rotating atomic gases[END_REF] and references therein). Due to the mathematical similarity between Coriolis and Lorentz forces, rotating neutral gases are indeed the exact analogue of an assembly of charged particles plunged in a magnetic field. Observing these highly correlated phases is one of the major goals in the field of trapped quantum gases [START_REF] Cooper | Rapidly rotating atomic gases[END_REF][START_REF] Bloch | Many-body physics with ultracold gases[END_REF]. This goal has, however, not yet been reached, because of the difficulty of communicating the required amount of angular momentum (on the order of N per atom, with N the number of particles) to the system [START_REF] Schweikhard | Rapidly rotating Bose-Einstein condensates in and near the lowest Landau level[END_REF][START_REF] Bretin | Fast rotation of a Bose-Einstein condensate[END_REF]. In practice, the residual static trap anisotropy limits the total angular momentum to much smaller values, for which the rotating gas is well described by a mean field approach [START_REF] Ho | Bose-Einstein condensates with large number of vortices[END_REF].

Recently alternative schemes to "simulate" artificial gauge fields for neutral atoms have been explored using two-dimensional (2D) optical lattices [START_REF] Jaksch | Creation of effective magnetic fields in optical lattices: the Hofstadter butterfly for cold neutral atoms[END_REF][START_REF] Ruostekoski | Particle number fractionalization of an atomic Fermi-Dirac gas in an optical lattice[END_REF][START_REF] Mueller | Artificial electromagnetism for neutral atoms: Escher staircase and Laughlin liquids[END_REF][START_REF] Sorensen | Fractional quantum Hall states of atoms in optical lattices[END_REF][START_REF] Dudarev | Spin-orbit coupling and berry phase with ultracold atoms in 2d optical lattices[END_REF][START_REF] Lim | Staggered-vortex superfluid of ultracold Bosons in an optical lattice[END_REF][START_REF] Klein | Phonon-induced artificial magnetic fields in optical lattices[END_REF]. As they do not involve any mechanical rotation of the system, they should be less sensitive to the imperfection of the trapping potential and thus easier to implement. The guideline for these proposals is the celebrated Harper model [START_REF] Harper | Single band motion of conduction electrons in a uniform magnetic field[END_REF][START_REF] Hofstadter | Energy levels and wave functions of Bloch electrons in rational and irrational magnetic fields[END_REF], defined by the two-dimensional (2D) single-particle Hamiltonian,

H Harper = -J n,m,± e ±i2παm ĉ †
n±1,m ĉn,m + ĉ † n,m±1 ĉn,m + h.c.

This tight-binding model was initially introduced [START_REF] Harper | Single band motion of conduction electrons in a uniform magnetic field[END_REF] to describe electrons hopping on a square lattice perpendicular to a constant magnetic field B. The operator ĉ † n,m creates a particle at position (x = nd, y = md), d is the lattice spacing, and J is the tunnel energy to nearest neighbors in the absence of magnetic field. The effect of a vector potential (here in the Landau gauge, A = -Bye x ) is included in the so-called Peierls phase factor, 2παm, where α = eBd 2 /h is the flux per unit cell expressed in units of the flux quantum φ 0 = h/e. The Harper Hamiltonian in Eq. ( 1) has been studied extensively in the literature (see, e.g. [START_REF] Harper | Single band motion of conduction electrons in a uniform magnetic field[END_REF][START_REF] Hofstadter | Energy levels and wave functions of Bloch electrons in rational and irrational magnetic fields[END_REF][START_REF] Thouless | Quantized Hall conductance in a two-dimensional periodic potential[END_REF][START_REF] Kohmoto | Zero modes and the quantized Hall conductance of the two-dimensional lattice in a magnetic field[END_REF]). For zero flux, the Harper model reduces to a standard tight-binding model on a square lattice, with a single Bloch band of width 8J. For rational fluxes α = p/q (where p and q are integers), the model is still periodic but with larger unit cells of size qd × d. The Bloch band splits into q magnetic sub-bands, and the ground state becomes q-fold degenerate. This results in the peculiar self-similar structure in an energy-magnetic flux diagram known as Hofstadter's butterfly [START_REF] Hofstadter | Energy levels and wave functions of Bloch electrons in rational and irrational magnetic fields[END_REF]. This structure exists in a strict sense for infinitely extended systems. For a finite system of size L, the details of this structure are washed out on fine scales such that q > L/d.

In their seminal proposition for an implementation of (1) [START_REF] Jaksch | Creation of effective magnetic fields in optical lattices: the Hofstadter butterfly for cold neutral atoms[END_REF], Jaksch and Zoller (JZ) suggested that the phase 2πα can be imprinted by a laser beam inducing hopping between adjacent sites (see also [START_REF] Ruostekoski | Particle number fractionalization of an atomic Fermi-Dirac gas in an optical lattice[END_REF][START_REF] Mueller | Artificial electromagnetism for neutral atoms: Escher staircase and Laughlin liquids[END_REF][START_REF] Sorensen | Fractional quantum Hall states of atoms in optical lattices[END_REF][START_REF] Dudarev | Spin-orbit coupling and berry phase with ultracold atoms in 2d optical lattices[END_REF][START_REF] Lim | Staggered-vortex superfluid of ultracold Bosons in an optical lattice[END_REF][START_REF] Klein | Phonon-induced artificial magnetic fields in optical lattices[END_REF] for alternative proposals). A major advantage of this approach is that large flux α ∼ 1 can be reached under experimentally realistic conditions. For a filling factor with ∼ 1 atom per site, the quantum Hall regime then becomes reachable in presence of moderate atomic interactions [START_REF] Sorensen | Fractional quantum Hall states of atoms in optical lattices[END_REF][START_REF] Palmer | High-field fractional quantum Hall effect in optical lattices[END_REF][START_REF] Palmer | Optical lattice quantum Hall effect[END_REF][START_REF] Hafezi | Fractional quantum Hall effect in optical lattices[END_REF][START_REF] Umucalilar | Trapped Fermi gases in rotating optical lattices: Realization and detection of the topological Hofstadter insulator[END_REF]. Theoretical studies have found that some of these phases were essentially similar to that encountered in the bulk [START_REF] Sorensen | Fractional quantum Hall states of atoms in optical lattices[END_REF][START_REF] Palmer | High-field fractional quantum Hall effect in optical lattices[END_REF][START_REF] Palmer | Optical lattice quantum Hall effect[END_REF][START_REF] Umucalilar | Trapped Fermi gases in rotating optical lattices: Realization and detection of the topological Hofstadter insulator[END_REF], whereas others are without counterpart in the bulk phase [START_REF] Hafezi | Characterization of topological states on a lattice with Chern number[END_REF][START_REF] Möller | Composite fermion theory for bosonic quantum Hall states on lattices[END_REF].

Although the existing proposals for implementing the Harper Hamiltonian with atomic quantum gases open many avenues, they are still challenging from an experimental point of view (for reasons discussed in detail below). In this paper, we extend the JZ proposal and discuss a new scheme based on an optical superlattice [START_REF] Sebby-Strabley | Lattice of double wells for manipulating pairs of cold atoms[END_REF][START_REF] Folling | Direct observation of second-order atom tunnelling[END_REF] to generate a gauge potential leading to [START_REF] Girvin | The quantum Hall effect: Novel excitations and broken symmetries[END_REF]. At variance with earlier works where alkali atoms were considered, we propose here to use atoms with a long-lived metastable excited state, such as alkaline-earth [START_REF] Kraft | Bose-Einstein condensation of alkaline earth atoms: 40 Ca[END_REF][START_REF] Stellmer | Bose-Einstein condensation of strontium[END_REF][START_REF] Martinez De Escobar | Bose-Einstein condensation of 84 Sr[END_REF] or Ytterbium [START_REF] Takasu | Spin-singlet Bose-Einstein condensation of two-electron atoms[END_REF][START_REF] Fukuhara | Degenerate Fermi gases of Ytterbium[END_REF] atoms. For concreteness we will discuss the case of Ytterbium, for which degenerate gases have already been produced for both bosonic and fermionic isotopes [START_REF] Takasu | Spin-singlet Bose-Einstein condensation of two-electron atoms[END_REF][START_REF] Fukuhara | Degenerate Fermi gases of Ytterbium[END_REF]. We take advantage of this level structure to alleviate many practical difficulties encountered The light-shift experienced by the atoms is α(λ)|E| 2 /2, with |E| the electric field strength at wavelength λ. Two particular wavelengths are of special interest, a "magic" wavelength λ m ≈ 760 nm and an "anti-magic" wavelength λ am ≈ 1.12 µm. (c) Sketch of the atomic configuration in a state-dependent optical lattice, showing the two interlaced sublattices for g (grey dots) and e (black dots) atoms. The y lattice is formed by a standing wave at the "magic" wavelength and confines both states identically. The x lattice is tuned near the "anti-magic" wavelength and confines the two states in two distinct sublattices shifted by λ am /4. with alkalis. Importantly, our scheme uses only building blocks which have already been individually demonstrated.

This paper is organized as follows. We first introduce the system, consisting of atoms with a ground state g and a long-lived metastable excited state e trapped in a twodimensional (2D) lattice. The optical lattice is tuned near an "anti-magic" wavelength where the polarisabilities of e and g are opposite [START_REF] Yi | State-dependent, addressable subwavelength lattices with cold atoms[END_REF]. The two internal states thus live in spatially separated sublattices. We show that the case α = 1/2, leading to Dirac points analogous to those appearing in graphene [START_REF] Castro Neto | The electronic properties of graphene[END_REF], can be achieved rather simply with a single laser beam coupling e and g. We then extend our scheme to generate arbitrary values of α using an optical superlattice, at the expense of several lasers operating at different frequencies. Subsequently, we address the influence of the additional potential that must be added to confine the particles. We investigate its effect by looking at (i) the density of states that determines the density profile of a trapped ideal Fermi gas, (ii) the ground state of the system, relevant for non-interacting bosons. In both cases, we obtain a clear signature for a non-zero value of α. We conclude by indicating how the method can be generalized to generate non-Abelian gauge fields [START_REF] Osterloh | Cold atoms in non-Abelian gauge potentials: From the Hofstadter "moth" to lattice gauge theory[END_REF].

Two-electron atoms in state-dependent optical lattices

Two-electron atoms, such as Ytterbium and alkaline-earth atoms, generally have a spinsinglet ( 1 S 0 ≡ g) ground state and an extremely long-lived spin triplet ( 3 P 0 ≡ e) excited state (see Fig. 1a). The very long lifetime of e (around 20 s for Yb [START_REF] Porsev | Possibility of an optical clock using the 6 1 S 0 → 6 3 P 0 transition in 171,173 Yb atoms held in an optical lattice[END_REF]) allows one to operate optical atomic clocks at the resonance wavelength λ 0 ≈ 578 nm between g and e [START_REF] Hong | Observation of the 1 S 0 -3 P 0 transition in atomic Ytterbium for optical clocks and qubit arrays[END_REF][START_REF] Hoyt | Observation and absolute frequency measurements of the S 1 -P 3 optical clock transition in neutral Ytterbium[END_REF][START_REF] Barber | Direct excitation of the forbidden clock transition in neutral 174 Yb atoms confined to an optical lattice[END_REF]. It has also inspired novel proposals for quantum information processing [START_REF] Reichenbach | Sideband cooling while preserving coherences in the nuclear spin state in group-II-like atoms[END_REF][START_REF] Stock | Entanglement of group-II-like atoms with fast measurement for quantum information processing[END_REF][START_REF] Gorshkov | Alkalineearth-metal atoms as few-qubit quantum registers[END_REF] or quantum simulation [START_REF] Sato | Non-Abelian topological order in s-wave superfluids of ultracold fermionic atoms[END_REF][START_REF] Gorshkov | Two-orbital SU(N) magnetism with ultracold alkaline-earth atoms[END_REF]. States g and e usually react differently to far-off resonance laser light, as they couple optically to many different states. The polarizabilities of the two states calculated using the data from [START_REF] Barber | Ytterbium Optical Lattice Clock[END_REF] are shown in Fig. 1b. For laser light tuned at the so-called magic wavelength λ m (∼ 760 nm for Yb), both states have the same polarizability and feel the same optical trapping potential [START_REF] Ye | Quantum state engineering and precision metrology using state-insensitive light traps[END_REF]. Central to our proposal is the existence of an "anti-magic" wavelength λ am ≈ 1.12 µm [START_REF] Yi | State-dependent, addressable subwavelength lattices with cold atoms[END_REF], for which the polarizabilities of g and e are opposite ‡ §. As we will see, this allows to create easily internal-state dependent potentials.

The first building block of our proposal consists in a laser beam at the magic wavelength, forming an horizontal light sheet that freezes the vertical motion of the atoms irrespective of their internal state e or g. From now on, we shall focus on the atom dynamics in the horizontal xy plane.The second ingredient consists in an optical standing wave along the y axis, also at the magic wavelength, which generates the same potential V (y) = -V 0 cos 2 (πy/d y ) for both states g and e (V 0 > 0). For the third building block, we apply another optical standing wave at the "anti-magic" wavelength λ am along the x axis, generating the potential V ′ (x) = ±V 1 cos 2 (πx/2d x ), where the + andsigns refer to e and g, respectively. The resulting potential V ′ (x) + V (y) creates two shifted sublattices (see Fig. 1c). Atoms in state g sit at positions r g = (2n, m) and atoms in state e at positions r e = (2n + 1, m), where n, m are integers. Here we take as lengths scales d x = λ am /4 along x, and d y = λ m /2 along y. We write the Wannier functions for the lowest Bloch band as w g (rr g ) and w e (rr e ). The potential V 0 is chosen such that the tunnelling energy J (y) along the vertical lines is significant. On the opposite, we choose V 1 large enough to neglect direct tunnelling along the x axis.

Tunnelling along the x axis is induced by a coherent optical coupling between g and e. Since atoms in e and g are trapped in separate optical lattices, this coupling induces hopping from one sublattice to the other. Following Jaksch and Zoller [START_REF] Jaksch | Creation of effective magnetic fields in optical lattices: the Hofstadter butterfly for cold neutral atoms[END_REF], we calculate the effective g -e hopping matrix element between two neighboring sites located at r g = (2n, m) and r e = r g + b = (2n + 1, m) as

J (x) eg e iq•rg = Ω 2 e iq•rg w * e (r -b)e iq•r w g (r) d 2 r, (2)
with Ω and q the laser Rabi frequency and wavevector. We suppose that the coupling laser propagates in the yz plane, so that q • r g = 2παm, where |α| can take any value between 0 and λ m /2λ 0 (≃ 0.66 for Yb) by adjusting the angle between the coupling laser beam and the z axis. ‡ There exists another anti-magic wavelength near 620 nm, which is unfortunately rather close to a transition connecting to the 3 P 0 state (649 nm). A lattice at this wavelength is thus likely to suffer from excessive spontaneous emission. § A many-electrons calculation [START_REF] Dzuba | Dynamic polarizabilities and related properties of clock states of Ytterbium atom[END_REF] using more accurate parameters for the positions and widths of the optical transitions predict the anti-magic wavelengths near 1120. Before proceeding with further analysis, let us give some numerical estimates for the different physical quantities involved in this experimental configuration. The anti-magic wavelength λ am corresponds for Yb to a recoil energy E R /h = h/(2mλ 2 am ) ≈ 900 Hz. Realistic values up to several tens of kHz can be achieved for the depth of the optical lattice potential V 1 /h. In order to avoid excitations to higher-lying Bloch bands, the Rabi frequency Ω must be small compared to ∆/ , where ∆ is the energy gap between the ground and first excited band. To give numerical estimates, we take V 1 = 20 E R , corresponding to a gap ∆ ≈ 8 E R . Assuming a moderate Rabi frequency Ω = E R / , we find a laser-induced tunneling energy J (x) eg ≈ 0.05 E R , much larger than the residual tunneling energy in each sublattice without change in internal state, J

(x) gg = J (x) ee ≈ 2.5 × 10 -3 E R .
In principle, besides resonant transitions between the fundamental bands of the two sublattices, the coupling laser also induces second-order tunneling processes, where an atom tunnels to a neighboring site in a given sublattice, via a virtual transition to an intermediate excited band in the other sublattice. We have verified that the corresponding tunneling rates were negligible compared to J (x) gg for the range of parameters investigated here. Finally, one can calculate the spontaneous emission and associated heating rates for atoms in each sublattice [START_REF] Gordon | Motion of atoms in a radiation trap[END_REF]. For lattice depths

V 0 = V 1 ∼ 30 E R , we find total heating rates γ (g) h ∼ 0.4 nK/s and γ (e) h ∼ 0.9 nK/s, respectively, or γ (g) h /E R ∼ 0.01 s -1 and γ (e)
h /E R ∼ 0.02 s -1 . The heating rates (dominated by the vertical arm at λ m ) are small enough to maintain the atoms in the ground band for several seconds. We first consider the simplest experimental arrangement, with a spin-dependent lattice as described above and a single coupling laser (Fig. 2a). The two-component g -e quantum gas, which is assumed to populate the lowest Bloch band only, is described by the Hamiltonian

H = -J n,m,n even e i2παm ĉ † 2n+1,m ĉ2n,m + e i2παm ĉ † 2n-1,m ĉ2n,m + h.c. -J n,m,± ĉ † n,m±1 ĉn,m + h.c. , (3)
when the lattice amplitudes are adjusted so that J (y) = J (x) eg = J. As in [START_REF] Jaksch | Creation of effective magnetic fields in optical lattices: the Hofstadter butterfly for cold neutral atoms[END_REF], the operator ĉ † n,m creates an atom at site (n, m), the parity of n identifying the internal state unambiguously. Alternatively, the internal degrees of freedom can be eliminated by going to a dressed state basis, and the spatial wavefunctions associated with each dressed state obey the same Harper equation [START_REF] Harper | Single band motion of conduction electrons in a uniform magnetic field[END_REF][START_REF] Hofstadter | Energy levels and wave functions of Bloch electrons in rational and irrational magnetic fields[END_REF] as the one derived from [START_REF] Girvin | The quantum Hall effect: Novel excitations and broken symmetries[END_REF]. Hence, at least in the absence of interactions (which we assume through this paper), the internal state is merely a label that can be ignored in the analysis.

In the original equation ( 1), the phase picked up across each link has the same sign for a given link direction in real space: e +i2παm in the +x direction and e -i2παm in the -x direction. Here, we find in contrast that the phase is tied to the direction in internal space: e +i2παm for g → e transitions and e -i2παm for e → g transitions. In real space, given the sublattices geometry, this corresponds to a phase that alternates in sign from an arbitrary unit cell (or "plaquette") to the neighboring one along the x axis. More precisely, consider a particle going clockwise around an elementary cell (figure 1c). The phase picked up along the path (2n, m -1) g → (2n, m) g → (2n + 1, m) e → (2n+1, m-1) e → (2n, m-1) g is the product of the two factors e i2παm e -i2πα(m-1) = e i2πα . Thanks to the phase factor ∝ e iq•rg , a particle making a loop around a plaquette acquires a non-zero global phase factor, reflecting the non-zero "magnetic flux" through the cell. However, for the adjacent plaquette, the path (also with clockwise orientation) (2n + 1, m -1) e → (2n + 1, m) e → (2n + 2, m) g → (2n + 2, m -1) g → (2n + 1, m -1) e leads to a phase factor e -i2πα opposite to the previous one. One then achieves in this way a staggered magnetic field, with a phase per plaquette ±2πα constant on a given column and changing sign from one column to the next (Fig. 2b).

The case α = 1 2 is peculiar, since in this case a staggered or a uniform flux both correspond to a tunneling phase e ±i2παm = (-1) m alternating between adjacent links depending on the parity of m. The Hamiltonian remains periodic, with a unit cell d x × 2d y and two sites per unit cell. The energy spectrum is given by

E(k x , k y ) = ±2J cos 2 (k x d x ) + cos 2 (k y d y ), (4)
where the Bloch vector (k x , k y ) sits in the first Brillouin zone ] -π/d x , π/d x ]×]π/2d y , π/2d y ]. We note that this dispersion relation gives rise to two 'Dirac points" for (k x , k y ) = (±π/2d x , π/2d y ) around which the dispersion relation is linear [START_REF] Lim | Staggered-vortex superfluid of ultracold Bosons in an optical lattice[END_REF][START_REF] Hou | Massless Dirac fermions in a square optical lattice[END_REF] (see also [START_REF] Zhao | BCS-BEC crossover on the two-dimensional honeycomb lattice[END_REF][START_REF] Zhu | Simulation and detection of Dirac fermions with cold atoms in an optical lattice[END_REF][START_REF] Wu | p x,y -orbital counterpart of graphene: Cold atoms in the honeycomb optical lattice[END_REF][START_REF] Lee | Ultracold fermions in a graphene-type optical lattice[END_REF][START_REF] Goldman | Non-Abelian optical lattices: Anomalous quantum Hall effect and Dirac fermions[END_REF] for discussions of the Dirac points occurring in different settings, such as hexagonal lattices or non-Abelian gauge fields). The physical consequences of the existence of these Dirac points are similar to those found in the case of Graphene [START_REF] Castro Neto | The electronic properties of graphene[END_REF]. However, their topological robustness [START_REF] Montambaux | Merging of Dirac points in a two-dimensional crystal[END_REF] might be affected by small deviations from α = 1/2 that break the two-site periodicity. For α = 1 2 , the staggered value of the phase per plaquette leads to a Hamiltonian that is different from (1), and which, in fact, does not exhibit any structure reminiscent of Hofstader's butterfly [START_REF] Wang | Hofstadter butterflies in macroscopic time-reversal invariant systems: Tight-binding electrons under n √ 2 × n √ 2 staggered magnetic fields[END_REF]. This can be linked to a general argument due to Haldane [START_REF] Haldane | Model for a quantum Hall effect without Landau levels: Condensed-matter realization of the "parity anomaly[END_REF], emphasizing that breaking time-reversal invariance symmetry is the key factor to generate systems with properties similar to (integer) quantum Hall phases. The situation considered up to this point, where the same laser beam drives all g → e transitions, is invariant with respect to time-reversal. We describe in the next section a mean to break this invariance and achieve a flux with a constant sign on all plaquettes.

State-dependent superlattice: Harper Hamiltonian for arbitrary α

Flux rectification

To ensure an identical phase 2πα per plaquette, instead of the staggered values ±2πα, we consider a situation where the on-site energies are modulated spatially along x with a period equal to twice the lattice spacing, as depicted in Fig. 2c (the y lattice remains identical for both states and is ignored in the following). Assume that the on-site energies are modulated according to

E g -V 1 , for x/d x = 4n, E e + V 2 , for x/d x = 4n + 1, E g -V 1 -V 2 , for x/d x = 4n + 2, E e , for x/d x = 4n + 3.
Here the E α (α = e, g) denote the internal energies in free space, where V 1 is the amplitude of the fundamental lattice with period d x and where V 2 ≪ V 1 is the amplitude of the modulation with period 2d x . In this potential landscape, the resonance frequencies for transitions linking neighboring sites become non-degenerate (Fig. 2c),

ω + = ω 0 + (V 1 + V 2 )/ for |g; 4n, m → |e; 4n + 1, m , ω -= ω 0 + (V 1 + 2V 2 )/ for |e; 4n + 1, m → |g; 4n + 2, m , ω + = ω 0 + (V 1 + V 2 )/ for |g; 4n + 2, m → |e; 4n + 3, m , ω ′ -= ω 0 + V 1 / for |e; 4n + 3, m → |g; 4n + 4, m ,
with ω 0 = E e -E g .

To the state-dependent lattice, one thus applies three coupling lasers propagating along y with frequencies ω + , ω -and ω ′ -. The laser at frequency ω + is chosen with a wavevector k L , and the lasers at frequencies ω -and ω ′ -with the opposite wavevector -k L . The alternance of the wavevectors compensates the alternance of the plaquette phases, thereby "rectifying" the staggered flux obtained in a single-frequency configuration. If we neglect off-resonant transitions, an analysis similar to that made in the previous section confirms that each plaquette now acquires the same phase factor e i2πα , and that the system is indeed described by the Harper Hamiltonian (1).

Optical superlattice

To realize such a modulation in practice, we propose to superimpose along the x direction an additional lattice potential with a double spatial period 2λ am , which can be written as V ′′ (x) = V g/e cos(πx/4d x +ϕ) 2 , where ϕ is the relative phase between the fundamental and double period lattices. Such superlattice potentials have been demonstrated experimentally for quantum gases [START_REF] Sebby-Strabley | Lattice of double wells for manipulating pairs of cold atoms[END_REF][START_REF] Folling | Direct observation of second-order atom tunnelling[END_REF]. This potential can be realized in practice by adding a second laser with frequency c/2λ am phase locked to the fundamental lattice to control the relative phase ϕ [START_REF] Folling | Direct observation of second-order atom tunnelling[END_REF]. The amplitudes V g/e depend on the internal state and are fixed by the initial choice of wavelength for the fundamental lattice. The analysis in Section 4.1 above (where V g = V e ) can be generalized straightforwardly to this less symmetric potential, with four non-degenerate transition frequencies in the general case. Moreover, it is possible to reduce the number of required laser frequencies to three by using a particular value of the relative phase ϕ that reproduces exactly the scheme described previously, with V 2 = V e cos ϕ -V g sin ϕ .

Conditions of validity

We now discuss the domain of validity of our model.

• First, as in the previous sections, we demand that laser-assisted tunneling dominates over regular tunneling without change in internal state. For a superlattice potential in the tight-binding regime, regular tunneling is suppressed when

V 2 ≥ E R , since J (x) gg , J (x)
ee ≪ E R . • Second, the ratio Ω/V 2 , with Ω the Rabi frequency characterizing the atom-laser coupling, must remain small to neglect off-resonant transitions and describe the Practically, the fluctuations δϕ of ϕ must be small enough to ensure that fluctuations of V 2 remain much smaller than the width of the ge resonance, set by the Rabi frequency Ω, or δϕ ≪ Ω/V 2 0.1. This is well within present technological capabilities.

system by an effective Harper model [START_REF] Jaksch | Creation of effective magnetic fields in optical lattices: the Hofstadter butterfly for cold neutral atoms[END_REF]. In practice, one can use a conservative value Ω = V 2 /10.

• Finally, the modulation V 2 must be small enough to avoid excitations to higher-lying Bloch bands by one of the coupling lasers. A reasonable choice is V 2 = ∆/3, where ∆ is the energy gap between the ground and first excited bands for V 2 = 0. This guarantees that transitions mixing different bands are detuned by at least ∆/3. For example, if we focus on the |e; 4n + 3, m → |g; 4n + 4, m transition, transitions linking the ground and first excited bands occur at frequencies ω ′ -± ∆/ . One can check that the detunings for the different coupling lasers are |V 2 ± ∆|, |2V 2 ± ∆|, ∆, which are indeed ≥ ∆/3 for the choice V 2 = ∆/3. We note that in the tight-binding regime, ∆ is notably larger than E R , which allows one to fulfill simultaneously all these conditions while maintaining a substantial value for J (x) eg . We summarize our considerations in Table 1, where we give typical values for lattice parameters and laser-assisted tunneling matrix elements in the superlattice configuration. We note that the values for J (x) eg remain sizeable, being on the same order as the tunnel energies J (x) gg in the g sublattice ( Ω = V 2 = 0) at moderate depths (∼ 10 E R ).

Depth Energy gap

Laser-assisted Table 1. Characteristic tunnel energies (along the x axis) achievable for realistic experimental parameters in the superlattice configuration. Here J eg denotes the laser-assisted coupling from one sublattice to the other. We used a Rabi frequency Ω = V 2 /10 = ∆/30 (see text).

V 1 /E R ∆/E R tunneling J (x) eg /E R

Comparison with the Jaksch-Zoller proposal

The scheme as described in this paper is inspired by the earlier proposal by Jaksch and Zoller [START_REF] Jaksch | Creation of effective magnetic fields in optical lattices: the Hofstadter butterfly for cold neutral atoms[END_REF]. That proposal was designed for alkali atoms in spin-dependent optical lattices. In this context, g and e denoted two different hyperfine states in the ground state manifold, and the spin-dependent optical potential was obtained by exploiting the vector light-shift arising in a laser field with suitable polarization (see [START_REF] Cohen-Tannoudji | Experimental study of Zeeman light shifts in weak magnetic fields[END_REF][START_REF] Mandel | Coherent transport of neutral atoms in spin-dependent optical lattice potentials[END_REF], for example).

The first difference between our method and the JZ proposal comes from the choice of two-electron atoms, versus "one-electron" alkalis. To obtain a useful spin-dependent lattice for akalis demands in practice to tune the lattice laser between the D 1 and D 2 lines, relatively close to resonance. For light alkalis, such as 6 Li or 40 K, the fine splitting between the D 1 and D 2 transitions is rather small (less than 1 nm), and spontaneous emission rates are too large to maintain the cloud in the nK regime. The most practical candidate seems to be 87 Rb, which has only stable bosonic isotopes and would probably still suffer from spurious spontaneous emission due to the coupling (Raman) lasers or from the lattice lasers. The choice of Yb removes all these drawbacks, by using only faroff resonant lasers and a dissipation-free coupling transition, and by allowing to study bosonic and fermionic systems in the same setup.

The use of a superlattice to ensure the same flux for all plaquettes constitutes the second main difference of our proposal with respect to JZ's proposal. There, a strong "tilting" potential independent of the internal state was introduced in order to distinguish |g; 2n, m → |e; 2n ± 1, m . Although this could be done in principle using a large static electric field or the light shift from an intense far-off resonant laser beam, it turns out to be rather challenging from an experimental point of view ¶. In addition, this tilting potential has to be almost perfectly linear over the extent of the cloud. Otherwise, the various transition frequencies would depend on space, spoiling the flux uniformity and the scalability of the method. We believe that the use of a superlattice optical potential will be a more flexible and more reliable tool that can be readily implemented experimentally.

Role of an external trap potential

We now discuss the properties of a system of non-interacting particles described by Eq. ( 1) in order to identify observable consequences of the laser-induced gauge potential. For a uniform system without trapping potential, and for α = 0, the model reduces to a standard tight-binding Hamiltonian. The spectrum consists in an allowed energy band of width 8J centered at zero energy. For the uniform Harper model, without trapping potential but with α = p/q = 0 (p and q are integers), the initial Bloch band splits into q magnetic sub-bands separated by well-defined energy gaps. This results in q distinct peaks in the density of state (DoS-see Fig. 3a andb). In a practical experiment, an additional trapping potential is inevitably present in addition to the optical lattice. Hence the physical system is both finite and inhomogeneous. This can strongly affect the properties of the uniform Harper model given in Eq. 1, but also provide novel experimental signatures, as we discuss in this section.

Density of states in presence of an external trap

We consider here the Harper model in Eq. ( 1) with an additional isotropic trapping potential V trap = κ(r/d) 2 . Let us recall that, with α = 0, one finds a quasi-continuum ¶ The Stark shift from the electric field is proportional to α|E| 2 , with α the static polarizability. Using a configuration where the cloud is offset from the center of the electrodes creating the field, one can create a linear potential of the form ηx/d x , with η ∼ αU 2 d x /D with D the distance between the electrodes and with U the applied electrostatic potential difference. For 87 Rb and D ∼ 1 cm, we find that a site-to-site energy offset on the order of the energy gap ∆, η ∼ ∆/3 ∼ 10 kHz requires U ∼ 50 kV, which is rather challenging to produce. of allowed states in the combined lattice plus external trap, in striking contrast with the uniform case [START_REF] Hooley | Single-atom density of states of an optical lattice[END_REF]. States with energy lower than 4J are qualitatively similar to the Bloch eigenstates of the uniform problem. These states are approximately contained in a disk of radius R eff ∼ J/κ. States with higher energy (which corresponds to the bandgap of the uniform system) are well-localized states offset from the cloud center [START_REF] Ott | Radio frequency selective addressing of localized atoms in a periodic potential[END_REF], with an energy mostly determined by the local potential energy near the turning point. In two dimensions, the density of states saturates at an asymptotic value of π/κ.

We now address the case with α = p/q = 0. The finiteness of the central region where tunneling is significant limits the relevant values of q to q < q max = R eff /d. We consider here the case α = 1/3 for concreteness. Solving the full model numerically, we find that the trapping potential has a dominant effect on the DoS and little remains of the initial peak structure (compare Fig. 3b and Fig. 3c). This indicates that global observables (such as thermodynamic quantities) will not be very sensitive to the presence of the gauge potential. We can understand this behavior from a semiclassical argument. Indeed, in a local density approximation, the density of states in a generic power-law potential κ(r/d) γ can be approximated as

ρ trap (E) = n d (2) rd (2) k (2π) 2 δ E -ǫ n (k) -κ r d γ (5)
where ǫ n (k) = 4Jχ n (k) is the relation dispersion for the nth magnetic sub-band of the untrapped system. This can be simplified into

ρ trap (E) = (4J) 1-γ 2 γκ n E 4J >χn(k) d (2) k S BZ E 4J -χ n (k) 2-γ . (6)
Here S BZ = (2π/d) 2 is the surface of the Brillouin zone. For a harmonic potential (γ = 2), the dependance on the reduced dispersion function χ is mild since it only arises from the upper bound of the integral: the peculiarities of the spectrum are essentially washed out by the spatial integration. This can be improved by using a steeper trap with γ > 2, e.g. a quartic potential as demonstrated in [START_REF] Bretin | Fast rotation of a Bose-Einstein condensate[END_REF].

Spatial distribution for non-interacting fermions

In the case of non-interacting fermions in a gauge potential, a more sensitive observable is given by the particle density [START_REF] Umucalilar | Trapped Fermi gases in rotating optical lattices: Realization and detection of the topological Hofstadter insulator[END_REF], assuming it can be measured in situ. Indeed, still within a semiclassical approximation, the spatial density can be written as

n(r) = n d (2) k (2π) 2 f FD [ǫ n (k) + V trap (r) -E F ] ,
with f FD the Fermi-Dirac distribution at T = 0 and E F the Fermi energy. We rewrite this as

n(r) = n d (2) k (2π) 2 ∞ -4J dE f FD [V trap + E -E F ] δ [ǫ n (k) -E] = ∞ -4J dE ′ f FD [V trap + E ′ -E F ] ρ hom [E ′ ] , (7)
where

ρ hom [E] = n d (2) k (2π) 2 δ [ǫ n (k) -E] (8)
is the density of state for the homogeneous system. For low enough temperatures, the clustering of energy levels in subbands in presence of a gauge potential will result in welldefined plateaux in the spatial density of the trapped cloud. These plateaux correspond to "band-insulating" regions, similar to that appearing for other incompressible systems in traps such as Mott insulators [START_REF] Folling | Formation of spatial shell structure in the superfluid to Mott insulator transition[END_REF][START_REF] Campbell | Imaging the Mott Insulator Shells by Using Atomic Clock Shifts[END_REF]. We have verified this fact in our numerical calculations, as shown in Fig. 3d. We conclude that the in-trap density profile will provide a clear signature of the gauge potential in the case of fermionic atoms.

Momentum distribution for non-interacting bosons

With non-interacting bosons, such a signature will be absent, as the atoms will condense into the lowest single-particle states (for the purpose of this paper, we do not consider interactions). However, signatures of the gauge potential will appear in reciprocal momentum space, which can be probed using standard time-of-flight imaging. The dephasing term e i2παm generates quasi-momentum components at the various harmonics The presence of a gauge potential reveals itself in momentum components at ±2π/3d in the former plot, and in the fine structures with period 3d in the latter. The confinement strength is κ/J = 0.005.

of the fundamental period 2πα. As a result, clear satellite peaks can be observed in the momentum distribution (see Fig. 4a), which will provide a direct experimental signature for bosonic species. We note in addition that the density distributions corresponding to these eigenstates shows periodic structures with period q (see Fig. 4b), although this signature appears on a much finer scale than in the fermionic case (thus being more difficult to observe).

Conclusion and perspectives

In this article, we have discussed how to implement Abelian artificial gauge potentials for Yb (or alkaline-earth) atoms, a scheme valid for bosons and fermions alike. In contrast to the bosonic isotopes with zero total spin, fermionic isotopes have nuclear spin. This potentially enables to simulate non-Abelian gauge fields, along the lines proposed by Osterloh et al. [START_REF] Osterloh | Cold atoms in non-Abelian gauge potentials: From the Hofstadter "moth" to lattice gauge theory[END_REF]. Let us focus for concreteness on the 171 Yb isotope, with nuclear spin I = 1 2 . In the presence of a moderate magnetic field of a few tens of Gauss, the degeneracy of the m I = ± 1 2 states within the e and g manifolds are lifted, giving in general four different internal transitions, ± 1 2 → ± 1 2 and ± 1 2 → ∓ 1 2 , with different frequencies. At the expense of using more laser frequencies and state-dependent lattices along both x and y directions, one can then arrange for the lasers to imprint different phases depending on the internal state or generate arbitrary rotations in internal space. As shown in [START_REF] Osterloh | Cold atoms in non-Abelian gauge potentials: From the Hofstadter "moth" to lattice gauge theory[END_REF], a straighforward generalization of the Abelian scheme then allows to engineer non-Abelian gauge potentials. Instead of a simple phase, hopping from sublattice g to e then corresponds to a rotation in internal space generated by noncommuting matrices. Combined with interatomic interactions, an entirely new class of superfluid or strongly correlated systems becomes accessible with ultracold atoms [START_REF] Osterloh | Cold atoms in non-Abelian gauge potentials: From the Hofstadter "moth" to lattice gauge theory[END_REF][START_REF] Goldman | Non-Abelian optical lattices: Anomalous quantum Hall effect and Dirac fermions[END_REF][START_REF] Goldman | Spatial patterns in optical lattices submitted to gauge potentials[END_REF][START_REF] Satija | Physics of a two-dimensional electron gas with cold atoms in non-Abelian gauge potentials[END_REF][START_REF] Goldman | Ultracold atomic gases in non-Abelian[END_REF]. In a wider context, this technique can be used to emulate spin-orbit coupling, similar to the Rashba interaction ∝ σx • k yσy • k x in planar semiconductors, or for the study of strongly interacting, ultracold fermionic matter in a non-Abelian background field, possibly relevant to high-energy physics.

We wish to thank the members of the "Bose-Einstein condensates" and "Fermi gases" group at LKB, G. Montambaux and G. Juzeliunas for helpful discussions, and A. Derevianko for a useful correspondence. This work was supported by ANR (Grant ANR-08-BLAN-65BOFL), IFRAF, the European Union (MIDAS project), and DARPA (OLE project). Laboratoire Kastler Brossel is unité mixte de recherche (UMR) n • 8552 of the CNRS.

Figure 1 .

 1 Figure 1. (a) Sketch of the level scheme showing the two internal states g and e coherently coupled by a laser beam; (b): Dynamic polarizability α(λ) of atoms in the ground (solid line) or excited states (dashed line).The light-shift experienced by the atoms is α(λ)|E| 2 /2, with |E| the electric field strength at wavelength λ. Two particular wavelengths are of special interest, a "magic" wavelength λ m ≈ 760 nm and an "anti-magic" wavelength λ am ≈ 1.12 µm. (c) Sketch of the atomic configuration in a state-dependent optical lattice, showing the two interlaced sublattices for g (grey dots) and e (black dots) atoms. The y lattice is formed by a standing wave at the "magic" wavelength and confines both states identically. The x lattice is tuned near the "anti-magic" wavelength and confines the two states in two distinct sublattices shifted by λ am /4.

3 .

 3 State-dependent lattice: Harper Hamiltonian for α = 1 2

Figure 2 .

 2 Figure 2. Laser schemes for realizing staggered (a,b) and uniform (c,d) magnetic fields. (a): The figure shows the different on-site energies along x, in the case of a simple-period lattice. All transition frequencies are degenerate, so that a single laser beam at frequency ω L is sufficient to couple all sites across the lattice. (b): Staggered effective magnetic field for a single coupling laser propagating along y with degenerate transition frequencies as shown in (a). (c): In a configuration where an additional, doubly-periodic shift of the on-site energies is introduced, the degeneracy between the transition frequencies is lifted. (d): With several dedicated laser beams (ω + propagating along y, ω -and ω ′-propagating along -y), one obtains a "rectified" uniform magnetic field.

Figure 3 .

 3 Figure 3. (a): Energy spectrum and (b): density of states (DoS) for the uniform model and α = 1/3. (c): density of states for the trapped model with α = 1/3 and confinement strength κ/J = 0.01. The dashed line shows the expected density of states for E ≫ J, irrespective of α. (d) Density profile for non-interacting fermions in a harmonic trap with κ/J = 0.01, for a flux α = 1/3 (solid lines) and Fermi energies E F = 0 (half-filled band) and 4J (filled band). The density profiles for α = 0 (dashed lines) are also shown for comparison. The calculation was performed on a 63×63 space grid.

Figure 4 .

 4 Figure 4. Typical momentum (a) and spatial (b) distributions in the ground state for non-interacting bosons and a flux α = 1/3.The presence of a gauge potential reveals itself in momentum components at ±2π/3d in the former plot, and in the fine structures with period 3d in the latter. The confinement strength is κ/J = 0.005.

gauge potentials: The case of constant Wilson loop. Phys. Rev. A, 79(2):023624, 2009.


Aude Ratier

Christelle Lopes

Gauthier Multari

Vanessa Mazerolles

Patrice Carpentier

Sandrine Charles
email: sandrine.charles@univ-lyon1.fr

New perspectives on the calculation of bioaccumulation metrics for active substances in living organisms

Keywords:

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

INTRODUCTION

Faced with the current environmental challenges, linked in particular to environmental pollution, ecotoxicology must today provide relevant and effective decision-making tools regarding bioaccumulation and effects of chemical substances on living organisms. Such tools must account for various exposure situations, environmentally realistic but complex (e.g., several routes of exposure, metabolism of substances, mixtures, etc.). Among available methods, toxicokinetic-toxicodynamic (TKTD) models are now strongly recommended to describe the link between exposure concentrations and effects on individual life-history traits over time from experimental data collected through toxicity tests, even standard ones (EFSA PPR Panel 2018). More specifically, the TK part of these models is used to relate the exposure concentration to the time course of the internal concentration within organisms, considering various processes such as accumulation, depuration, metabolization and excretion (known as ADME processes). As some recent regulations, the EU regulation No 283/2013 for plant protection products in marketing authorisation applications requires for example a bioaccumulation test on fish according to OECD Test guideline 305 (OECD 2012), which consists in an accumulation phase followed by a depuration phase. During the accumulation phase, fish are exposed to a substance of interest at a range of concentrations, chosen according to the assumed mode of action of the substance. After a certain time period fixed by the experimenter, organisms are transferred to a clean medium for a depuration phase. The concentration of the substance (and of its potential metabolites) is followed within fish at regular time points during both phases leading in fine to the estimation of bioaccumulation metrics. In this paper, for the sake of generality, we chose the generic expression "bioaccumulation metrics" to denote either bioconcentration factors (BCF) used when exposure is via water, biota-sediment accumulation factors (BSAF) when exposure is via sediment or biomagnification factors (BMF) when exposure is via food. Bioaccumulation metrics appeared to us as the best compromise regarding the wide diversity of terms used in the scientific literature (for example, USEPA (1994), [START_REF] Gobas | Revisiting Bioaccumulation Criteria for POPs and PBT Assessments[END_REF] and [START_REF] Burkhard | Comparing laboratory and field measured bioaccumulation endpoints[END_REF]).

Nevertheless, preference is often given to experimentally derived BCF estimates to be used for secondary poisoning assessment under Biocidal Products Regulation (European Commission 2012). Bioaccumulation tests are of course not only limited to fish, even if a test according to OECD Test guideline 305 (OECD 2012) is preferred when experimental information on bioaccumulation is needed for PBT/vPvB assessment under REACH regulation (ECHA 2017;European Commission 2006). Consequently, any other species can be used, such as benthic invertebrates, terrestrial oligochaetes or birds, depending on the substance under consideration.

From a regulatory point of view, bioaccumulation metrics are key decision criteria used to evaluate concentrations of active substances in food items of vertebrates (especially piscivorous birds and mammals), making the estimation of these metrics with the most precision as possible a highly crucial methodological challenge.

All bioaccumulation metrics rely on estimates of kinetic parameters as involved in toxicokinetic (TK) models. In the past decades, many types of methods have been proposed to get these estimates from simple TK models, most of them providing BCF estimates separately considering the kinetics for both accumulation and depuration phases as observed in dedicated experiments (OECD 2012). Nevertheless, TK model parameters are known to be highly correlated, so that separating their estimation prevents to account for a mutual influence on their uncertainty. Moreover, kinetic parameter estimates are only rarely provided with their uncertainty, although this is now expected by the regulatory bodies (EFSA Scientific Committee 2018). Consequently, concomitantly to the above-mentioned challenge, environmental risk assessment could be improved if complete tools allowing for a simultaneously estimation of all TK model parameters associated with their uncertainty would be available in support of stakeholders that need to fulfil regulatory expectations. These tools must also be easy to use to overcome the scepticism of regulators who are faced with multiple TK models and implementation methods, while thinking about their standardization at the same time (Tan et al. 2020). An R-code was first proposed in 2016 [START_REF] Aldenberg | bmcfR: Tools for Modeling Bioaccumulation Potential in Fish[END_REF]) that allowed to analyse data collected only from the OCED test guideline 305. In the same line of thought, a spreadsheet was recently proposed by [START_REF] Gobas | A Toxicokinetic Framework and Analysis Tool for Interpreting Organisation for Economic Co-operation and Development Guideline 305 Dietary Bioaccumulation Tests[END_REF]. Our paper goes beyond by considering all type of species-compound accumulation-depuration data (not only those of the OECD test guideline 305) which analysis leads to one or several bioaccumulation metrics of interest. [START_REF] Ratier | A Bayesian framework for estimating parameters of a generic toxicokinetic model for the bioaccumulation of organic chemicals by benthic invertebrates: Proof of concept with PCB153 and two freshwater species[END_REF] recently proposed a full revisit of the TK modelling approach based on a unified inference method to estimate parameters of TK models for both accumulation and depuration phases, simultaneously, automatically associating the uncertainties. This innovative framework has been thought to make it possible to further incorporate the TK part into complete TKTD models. Benefiting from this innovation, we present new perspectives for a facilitated calculation of any type of bioaccumulation metrics (such as BCF/BSAF/BMF) thanks to the new ready-to-use statistical web tool MOSAICbioacc (http://mosaic.univ-lyon1.fr/bioacc).

MOSAICbioacc runs one-compartment TK models that are automatically designed according to the input data. MOSAICbioacc leads to bioaccumulation metrics associated with their uncertainty propagated from the kinetic parameter estimates, without the need for users to invest underlying technical aspects. MOSAICbioacc is free of use, fully integrated within the all-in-one facility MOSAIC itself (http://mosaic.univ-lyon1.fr). MOSAICbioacc is regularly updated to always offers the very latest conceptual advances related to TK models. Today, MOSAICbioacc allows accounting for several exposure routes (water, pore water, sediment and food), for metabolism of chemicals (if the input experimental data include measurements for both the parent chemical and its metabolites), and for potential growth of organisms (if growth measurements are included within the data set). The use of MOSAICbioacc only requires users to upload their experimental data, collected via standard protocols or from home-made experimental designs.

Priority was first given to the calculation of BCF/BSAF/BMF because they are the widely used bioaccumulation metrics in the current regulatory guidelines. MOSAICbioacc get their estimate as probability distributions that are summarized for users by the median (50 th centile of the distribution) and the 95% credible interval (delimited by the 2.5 th and 97.5 th centiles) quantifying the uncertainty. In addition, the fitting plots and all model parameter estimates are provided, followed by a collection of goodness-of-fit criteria allowing users to check the relevance of the results. All outputs can be downloaded under different formats for further inclusion into any home-made document (in particular the open-source programming code), and a full report can also be directly downloaded gathering everything that is displayed to users on the web page. These two latter features of MOSAICbioacc guarantee both reproducibility and transparency of underlying calculations.

MATERIALS AND METHODS

MOSAICbioacc is part of the web platform MOSAIC (https://mosaic.univ-lyon1.fr/, [START_REF] Charles | MOSAIC: a web-interface for statistical analyses in ecotoxicology[END_REF]. It was developed as a Shiny environment [START_REF] Cheng | Shiny: Web Application Framework for R[END_REF]), available at https://mosaic.univ-lyon1.fr/bioacc and hosted at the Rhône-Alpes Bioinformatics Center PRABI (PRABI 2020). A user guide and an explanatory video are immediately available in the introductory section of the application to fully assist users step by step in appropriating the tool and its features. Details on underlying ordinary differential equations and their solving are also provided within a detailed user guide (see supplemental data, Annex 1), thus ensuring all required transparency as recommended by EFSA for a good modelling practice (EFSA PPR Panel 2014).

Data uploading

When using MOSAICbioacc, the first step is to upload input data (Fig. 1-a). MOSAICbioacc expects to receive experimental exposure time-course data, including at least an accumulation phase, as a .txt file or a .csv file (comma, semicolon or tabular separator) with a specific format.

Each line of the table corresponds to one time point for a given replicate and a given exposure concentration of the contaminant. The data table must contain at least four columns, with the exact following headers, the order of columns being not important (Table 1): 'time' (the time point of the measurement at the exposure concentration, in hours, minutes, days or weeks); 'expw', 'exppw', 'exps' or 'expf': the exposure concentration of the contaminant in the medium, that is water, pore water, sediment or food, respectively, all expressed in µg.mL -1 or in µg.g -1 ); 'replicate' (a number or a character that is unique for each replicate, dimensionless); and 'conc' (the concentration of the contaminant, and of its potential metabolites, measured within organisms, in µg.g -1 ). According to the experimental design, further columns can be added in the data file: 'expw', 'exppw', 'exps' and/or 'expf' if several exposure routes are considered together, 'concmℓ' (the concentrations of metabolite ℓ derived from the parent compound within the organisms, e.g., concm1, concm2, ...), and 'growth' (if growth measurements of the organisms are available).

As shown on Fig. 1-a, MOSAICbioacc users can either upload their own data set with a click on 'Browse' (by taking care about the expected format specification) or try MOSAICbioacc with example data sets (in total six data sets are proposed, each with different characteristics). When the upload is complete, users must manually select the appropriate separator, the time unit and the duration of the accumulation phase; please note that, when using example data sets, these fields are automatically filled in.

MOSAICbioacc first provides a table with the raw data allowing users to check if the data were correctly entered. Users can also visualize a plot with the raw data (Fig. 1-b). In case data were collected for several exposure concentrations and if the users have uploaded all, one must be chosen for the MOSAICbioacc analysis. Note that only one file at a time can be analysed. Also, when another exposure concentration from the same data file is chosen, the duration of the accumulation phase is reset and framed in orange to invite users for update before to launch new calculations. Example files provided in MOSAICbioacc are dedicated to assist users in formatting their own data and to appropriate the different MOSAICbioacc features from the very first step. This paper illustrates this step-by-step process based on a typical data set of a toxicokinetic bioassay where internal concentrations were collected in fathead minnows (Pimephales promelas) exposed via contaminated water to a highly hydrophobic chemical (logKow = 9.06) at an exposure concentration of 0.0044 µg.mL -1 over 49 days, with one replicate at each time point. After 49 days, minnows were transferred in a clean medium for 98 days more [START_REF] Crookes | Estimation of fish bioconcentration factor (BCF) from depuration data[END_REF]. The data set (CSV format) and the report with all results (HTML format) can be downloaded directly from MOSAICbioacc.

Model and parameters

All TK models considered in MOSAICbioacc describe organisms as single compartments for which a first-order kinetic bioaccumulation model accounting for several exposure routes and elimination processes can be expressed in a generic way as follows (Eqs. ( 1) to ( 4)):

⎩ ⎨ ⎧ 𝑑𝐶 (𝑡) 𝑑𝑡 = 𝑈 -(𝐸 + 𝑀)𝐶 (𝑡) (1) 𝑑𝐶 ℓ (𝑡) 𝑑𝑡 = 𝑘 ℓ 𝐶 (𝑡) -𝑘 ℓ 𝐶 ℓ (𝑡) (2) for 0 ≤ 𝑡 ≤ 𝑡 ⎩ ⎨ ⎧ 𝑑𝐶 (𝑡) 𝑑𝑡 = -(𝐸 + 𝑀)𝐶 (𝑡) (3) 𝑑𝐶 ℓ (𝑡) 𝑑𝑡 = 𝑘 ℓ 𝐶 (𝑡) -𝑘 ℓ 𝐶 ℓ (𝑡) (4)
for 𝑡 > 𝑡

where Cp(t) is the internal concentration of the parent compound at time t (µg.g -1 ), 𝐶 ℓ (t) the internal concentration (µg.g -1 ) of metabolites (∀ℓ = 1 … 𝐿 with L the total number of metabolites) at time t, U the sum of all uptake terms, E the sum of all elimination terms for the parent compound, M the sum of all metabolization terms, 𝑘 ℓ the metabolization rate of metabolite ℓ (time -1 ) and 𝑘 ℓ the elimination rate of metabolite ℓ (time -1 ). Table 2 gives an overview of all parameter and variable meaning. The dynamical system in equations ( 1) to (4), corresponding to the deterministic part of the model, can explicitly be solved when the exposure concentration is assumed to be constant over time (Annex 1). A Gaussian probability distribution was assumed as the stochastic part of the final model, based on the quantitative continuous nature of the concentration variables for both the parent compound and its metabolites within the organisms (Eqs. ( 5) and ( 6)):

𝐶 , (𝑡) ∼ 𝒩 𝐶 (𝑡), 𝜎 (5)
𝐶 , ℓ (𝑡) ∼ 𝒩 𝐶 ℓ (𝑡), 𝜎 ℓ (6)

Today MOSAICbioacc proposes data analyses by including until four exposure routes via water, pore water, sediment, and/or food), until three processes of elimination, which are excretion, growth dilution and biotransformation, with a maximum of 15 metabolites directly deriving from the parent compound (i.e., phase I metabolism). For example, if three metabolites are considered (L = 3), the number of parameters involved in the most complete TK model equals 19. In total, users have 112 possible models, automatically designed according to their data (see supplemental data, Annex 2). For a given data set, the most complete model is built up by default and first proposed to users, from which they can perform the MOSAICbioacc analysis (Fig. 1-c). Users can also deselect some of the parameters (based on biological hypotheses related to the most probable exposure route or by neglecting one elimination process, for example). These choices lead to the automatic building of a nested TK sub-model to fit again on the data.

For clarity reasons, this paper illustrates MOSAICbioacc features from a simple data set (Pimephales_two.csv file) considering only the water exposure route (parameter kuw for the uptake rate) and the excretion process (parameter kee for the elimination rate), that is one of the simplest TK models among the 112 possibilities. The corresponding equations for the deterministic part of this TK model are given below (Eqs. ( 7) and ( 8)):

⎩ ⎨ ⎧ 𝑑𝐶 (𝑡) 𝑑𝑡 = 𝑘 × 𝑐 -𝑘 × 𝐶 (𝑡) for 0 ≤ 𝑡 ≤ 𝑡 (7) 𝑑𝐶 (𝑡) 𝑑𝑡 = -𝑘 × 𝐶 (𝑡) for 𝑡 > 𝑡 (8)
Finally, there are only three parameters to estimate: kuw, kee and parameter σCp of the Gaussian distribution (Eq. ( 5)). An illustration of a more complex data set with biotransformation and growth processes for a benthic invertebrate is given as supplemental material (Annex 3), as well as examples for exposure route by sediment (Annex 4) or food (Annex 5).

Once the Pimephales_two.csv example file has been uploaded, all required fields related to the experimental design are automatically filled in and raw data can be visualized, either as a table (default) or a plot, and the button to launch calculations ('Calculate and Display') is unlocked (Fig. 1-c). Once this button clicked, calculations start running with a progress bar informing users about progress. When calculations are finished, results are displayed, as plots or tables first for the bioaccumulation metrics, then for the fitting plots and some relevant goodness-of-fit criteria.

Calculations

Bayesian inference

Computations underlying MOSAICbioacc results are performed with JAGS [START_REF] Plummer | rjags: Bayesian Graphical Models using MCMC[END_REF]) and the R software (R Core Team 2020, version 4.0.2) via the rjags and jagsUI packages [START_REF] Plummer | rjags: Bayesian Graphical Models using MCMC[END_REF][START_REF] Kellner | jagsUI: A Wrapper Around 'rjags' to Streamline[END_REF]. Models are fitted to bioaccumulation data using Bayesian inference via Monte Carlo Markov Chain (MCMC) sampling. For each model, calculation running starts with a short sampling on three MCMC chains (5,000 iterations after a burn-in phase of 10,000 iterations) using the Raftery and Lewis method [START_REF] Raftery | Practical Markov chain Monte Carlo]: Comment: one long run with diagnostics: implementation strategies for Markov chain Monte Carlo[END_REF] to set the necessary thinning and the appropriate number of iterations in order to reach a precise and accurate estimation of each model parameter. Thanks to rjags, model parameters are retrieved as a joint posterior distribution from the likelihood of the observed data combined with prior information on parameters. All details on this approach can be found in the original research paper [START_REF] Ratier | A Bayesian framework for estimating parameters of a generic toxicokinetic model for the bioaccumulation of organic chemicals by benthic invertebrates: Proof of concept with PCB153 and two freshwater species[END_REF]) but also in many other papers in the field of ecotoxicology [START_REF] Billoir | Bayesian modelling of daphnid responses to time-varying cadmium exposure in laboratory aquatic microcosms[END_REF]EFSA PPR Panel 2018).

Choice of prior distributions

For simplicity reasons, we hid the choice of priors to MOSAICbioacc users; hence, they cannot be changed, except by downloading the open source programming code and handling it directly within the R software. To ensure genericity of priors we chose non-informative (-5, 5) log10-uniform distributions for all uptake and elimination rate constants, and non-informative (0, A) uniform distributions for all standard deviations with a large A, here defined as five times the maximum internal measured concentration, which is then removed from the data set, as usually proceeded [START_REF] Gelman | Prior Distribution for Variance Parameters in Hierarchical Models[END_REF]).

Bioaccumulation metrics

Bioaccumulation metrics are the first outputs delivered by MOSAICbioacc. From the example chosen for this paper, the data analysis led to both the kinetic bioconcentration factor (BCFk) and the steady state bioconcentration factor (BCFss), with the following exact mathematical expressions (Eqs. ( 9) and ( 10)):

𝐵𝐶𝐹 = 𝑘 𝑘 (9) 𝐵𝐶𝐹 = 𝐶 (𝑡 ) 𝑐 (10)
where 𝐶 (𝑡 ) is the internal parent compound concentration (in µg.g -1 ) at the end of the accumulation phase (that is at t = tc, in time) and 𝑐 is the exposure contaminant concentration in water (µg.mL -1 ). More details about calculations of bioaccumulation metrics at steady state are provided in Annex 6.

RESULTS AND DISCUSSION

Bioaccumulation metrics

When users click on the 'Calculate and Display' button, some results are provided by default. First, the BCFk is given as a probability distribution (Fig. 2-a) and summarized with its median and its 95% uncertainty limits, that is 95% credible interval delimited by the 2.5 th and the 97.5 th centiles of the posterior probability distribution (Table 3-a). If users ask for the BCFss, its probability distribution is also delivered (Fig. 2-b) also summarized with the median and the 95% uncertainty limits (Table 3 -a). Credible intervals are crucial information to quantify the uncertainty on parameter estimates. If data are available for several exposure routes (according to the experimental design) and uploaded within MOSAICbioacc, the BCF/BSAF/BMF metrics are displayed in separate tabs (see Annex 4 for an example).

Predictions

Following bioaccumulation metrics calculations, the fitted curve and its uncertainty band superimposed to the observations is provided (Fig. 2-c). From the joint posterior distribution of model parameters, MOSAICbioacc then provides the marginal posterior distributions for each parameter, which are also summarized with quantiles in a table (Table 3- b): medians (for point estimates) and 2.5 th and 97.5 th centiles (for 95% credible intervals).

Goodness-of-fit criteria

After fitting plots, several goodness-of-fit criteria follow in a prioritized order chosen based on their relevance and their ease of interpretation. The fitting quality of the model can be first checked using the Posterior Predictive Check (PPC) plot: the idea is to compare each observed value to its prediction from the fitted model at the corresponding exposure concentration associated with its 95% credible interval. If the fit is correct, we expect to get 95% of the observed values falling within the 95% credible intervals of their predictions. As shown on Fig. 3-a, x-axis locates the observed values, while the y-axis reports their median predictions (black dots) with their 95% credible intervals (vertical segments).

The relevance of the inference process can also be checked using the comparison of prior and posterior distributions for each model parameter. The overall expectation is to get a narrower posterior distribution compared to the prior one for each parameter, reflecting that data contributed enough to precisely estimate parameters (Fig. 3-b). Users have the possibility to select plots for deterministic (e.g., kuw, kee) or stochastic (e.g., σCp) parameters.

Then, MOSAICbioacc provides a coloured matrix in order to see at a glance the most correlated or anti-correlated parameters, in order to quickly diagnose potential problems of precision due to highly correlated parameters. Moreover, MOSAICbioacc provides plots to visualize correlations between parameters (Fig. S1-a). Such a plot is obtained by projecting the joint posterior distribution as a matrix in planes of parameter pairs where contours have shapes reflecting both the sign and the strength of the correlations (sub-diagonal). The correlation plot also gives marginal posterior distribution of each model parameter (diagonal) and Pearson correlation coefficients (upper diagonal). Correlations between parameters are important to consider in particular when they are high (namely, greater than 0.75) what would mean that one parameter estimate could considerably influence the other, and reciprocally. Users can display the correlation plot for deterministic parameters only or for all parameters.

The convergence of MCMC chains can be checked with the Gelman-Rubin diagnostic [START_REF] Gelman | Inference from Iterative Simulation Using Multiple Sequences[END_REF] expressed via the potential scale reduction factor (PSRF) which is expected to be close to 1.00 (Fig. S1-b). It can be also visually verified from the MCMC trace plots, which show the time series of the sampling process leading to the posterior distribution for each parameter; it is expected to get overlap of all MCMC chains (Fig. S1-c). Users can visualize the MCMC trace plots for deterministic or stochastic parameters.

Finally, the Deviance Information Criterion (DIC) is provided. It is a penalized deviance statistics accounting for the number of parameters that is only useful to compare several models fitted to a same data set. Models with lower DIC values will be preferred. So, the DIC is only useful when several sub-models are compared based on different choices of parameters from all the possible combinations that the users can choose from the beginning according to the uploaded data.

Downloads

At the bottom of the result web page, all outputs, either separately or as a full report, can be downloaded. Users can also download the entire R code corresponding to all calculations and graphs from the uploaded data set. This ensures transparency and reproducibility of all MOSAICbioacc results. This R code can be used as a steppingstone to change default options or to perform further analyses directly in the R software. For example, users can modify figures at their convenience or make several analyses on several data sets at the same time.

CONCLUSION

Offering MOSAICbioacc as a new on-line service, free, user-friendly and ready-to-use, raises important methodological issues: (1) automation of the inference process, in particular for Bayesian inference;

(2) options and choices in a transparent and facilitated way for users;

(3) default outputs whose order has been chosen to facilitate their step-by-step interpretation;

(4) easy accessibility to figures, tables and R code to be downloaded under different convenient formats; and (5) a final full report of all MOSAICbioacc analyses. Besides its user-friendliness, MOSAICbioacc is free of use while ensuring privacy of the uploaded data as well as transparency and reproducibility of results, together with a short response time. MOSAICbioacc is particularly useful to estimate parameters of TK models leading to predictions of chemical concentrations bioaccumulated within living organism (whatever the species, aquatic, aerial or terrestrial) from accumulation-depuration data, even standard ones. MOSAICbioacc could thus be of particular interest for risk assessors and decision makers in their daily work of evaluating dossiers, e.g., for market authorisation of active substances (EFSA PPR Panel 2018). Indeed, all results provided by MOSAICbioacc account for uncertainty and correlations between parameters, making possible to reproduce any previous analysis that would need to be confirmed.

MOSAICbioacc can also be used for more exploratory research purposes by any environmental scientists or ecotoxicologists when accumulation-depuration data are collected and need to be Observed (black dots) and predicted contaminant concentrations in the organisms (µg.g -1 ),

where the median curve is displayed as the solid orange line and the uncertainty band as the grey zone, delimited by the 2.5 th and 97.5 th centiles in dotted orange lines.

LIST OF TABLES

Table 1. Example of a data set ready to be used in MOSAICbioacc. The data set must contain four columns whatever their order: 'time' (time points of the measurements at the exposure concentration, in time unit: here in days), 'conc' (contaminant concentrations measured within organisms, that must be expressed in µg.g -1 ), 'expw' (the contaminant concentration in the exposure medium, here water, that must be expressed in µg.mL -1 ), and 'replicate' (a number or a character that is unique for each replicate, dimensionless). Table 1. Example of a data set ready to be used in MOSAICbioacc. The data set must contain four columns whatever their order: 'time' (time points of the measurements at the exposure concentration, in time unit: here in days), 'conc' (contaminant concentrations measured within organisms, that must be expressed in µg.g -1 ), 'expw' (the contaminant concentration in the exposure medium, here water, that must be expressed in µg.mL -1 ), and 'replicate' (a number or a character that is unique for each replicate, dimensionless).

Figure 1 .

 1 Figure 1. Data uploading and user information as required from the MOSAICbioacc homepage.

Figure 2 .

 2 Figure 2. Probability distributions of BCFk (a) and BCFss (b). The middle-dotted line represents

Figure 3 .

 3 Figure 3. Some goodness-of-fit criteria as provided by MOSAICbioacc: (a) Posterior Predictive

Figure 1 .Figure 2 .Figure 3 .

 123 Figure 1. Data uploading and user information as required from the MOSAICbioacc homepage.

  (𝑡) internal measured concentration of the parent compound at time 𝑡 (in 𝜇𝑔. 𝑔 -1 ) 𝐶 𝑜𝑏𝑠,𝑚 ℓ (𝑡)internal measured concentration of metabolite ℓ (in 𝜇𝑔. 𝑔 -1 ) 𝜎Cp standard deviation for internal concentration of the parent compound (in 𝜇𝑔. 𝑔 -1 ) 𝜎 𝑚𝑒𝑡 ℓ standard deviation for the internal concentration of metabolite ℓ (in 𝜇𝑔. 𝑔 -1 ) elimination terms for metabolite ℓ

Table 2 .

 2 Meaning of parameters and variables of the TK model used in MOSAICbioacc.

Table 3 .

 3 Example of BCF (a) and parameter (b) estimates expressed as medians (50 th centile) and 95% credibility intervals (2.5 th -97.5 th centiles). Hyphens stand for dimensionless parameters.Tan Y.M., Chan M., Chukwudebe A., Domoradzki J., Fisher J., Hack C.E., Hinderliter P., Hirasawa K., Leonard J., Lumen A., Paini A., Qian H., Ruiz P., Wambaugh J., Zhang F., Embry M. 2020. PBPK model reporting template for chemical risk assessment applications. Office of water. Off. of Sci. and Tech., Washington, DC, USA.

	Regul. Toxicol. Pharmacol. 115. doi:10.1016/j.yrtph.2020.104691
	USEPA. [USEPA] U.S. Environmental Protection Agency. 1994. Great Lake water
	quality initiative technical support document for the procedure to determine bioaccumulation
	factors 822-R-94-002.

Table 2 .

 2 Meaning of parameters and variables of the TK model used in MOSAICbioacc. 𝑐 𝑖 exposure concentration of route 𝑖 (in 𝜇𝑔. 𝑚𝐿 -1 ) 𝐶 𝑝 (𝑡) internal concentration of the parent compound at time 𝑡 (in 𝜇𝑔. 𝑔 -1 ) 𝐶 𝑚 ℓ (𝑡) internal concentration of metabolite ℓ (in 𝜇𝑔. 𝑔 -1 ) 𝑘 𝑢 𝑖 uptake rate of exposure route 𝑖 (expressed in 𝑡𝑖𝑚𝑒 -1 ) 𝑘 𝑒 𝑗 elimination rate of elimination process 𝑗 (expressed in 𝑡𝑖𝑚𝑒 -1 ) 𝑘 𝑒 𝑚ℓ elimination rate of metabolite ℓ (expressed in 𝑡𝑖𝑚𝑒 -1 ) 𝑘 𝑚 ℓ biotransformation rate of metabolite ℓ (expressed in 𝑡𝑖𝑚𝑒 -1 ) 𝑡 𝑐 duration of the accumulation phase (in time unit) 𝒩 Gaussian probability distribution 𝐶 𝑜𝑏𝑠,𝑝

	time	conc	expw	replicate
	0	0.000	0.0044	2
	3	0.225	0.0044	2
	7	0.355	0.0044	2
	14	0.553	0.0044	2
	21	0.658	0.0044	2
	28	0.785	0.0044	2

Table 3 .

 3 Example of BCF (a) and parameter (b) estimates expressed as medians (50 th centile) and 95% credibility intervals (2.5 th -97.5 th centiles). Hyphens stand for dimensionless parameters.

		2.5 th	50 th	97.5 th	Units
	(a)
	BCFk	233	276	338	-
	BCFss	147	233	317	-
	(b)
	kuw	7.388	10.59	15.45	d -1
	kee	0.02307	0.03851	0.06091	d -1
	σCp	0.1249	0.1684	0.2445	µg.g -1

The authors are thankful to ANSES for providing the financial support. The MOSAICbioacc web tool is hosted at the Rhône-Alpes Bioinformatics Center PRABI (PRABI, 2020). This work benefited from the French GDR "Aquatic Ecotoxicology" framework which aims at fostering stimulating scientific discussions and collaborations for more integrative approaches. This work is part of the ANR project APPROve (ANR-18-CE34-0013) for an integrated approach to propose proteomics for biomonitoring: accumulation, fate and multimarkers (https://anr.fr/Projet-ANR-18-CE34-0013). This work was also made under the umbrella of the Graduate School H2O'Lyon (ANR-17-EURE-0018) and "Université de Lyon" (UdL), as part of the program "Investissements d'Avenir" run by "Agence Nationale de la Recherche" (ANR). The authors are truly grateful to the anonymous colleagues who participated in testing MOSAICbioacc and for giving their feedback. The authors also are grateful for Benoît BRET for creating the logo for MOSAICbioacc. The authors declare no

accessible directly within MOSAICbioacc at https://mosaic.univ-lyon1.fr/bioacc.

available as a new statistical service dedicated to TK modelling approaches embedded within the MOSAIC platform [START_REF] Charles | MOSAIC: a web-interface for statistical analyses in ecotoxicology[END_REF]. This makes MOSAIC an all-in-one facility for many applications. MOSAICbioacc will be shortly extended with a prediction tool, in order to help in designing new TK experiments, for example, when a new species-compound combination requires attention or when additional data are needed to get a better precision on bioaccumulation metrics. To gain again in generality, an R-package (rbioacc) is already in progress to include all the functionalities encompassed in MOSAICbioacc, which will also make it possible to extend the use of TK models when the exposure concentrations vary over time, and when it is necessary to consider non-first order kinetics.

SUPPLEMENTAL DATA

Supplemental data are available on-line.