Chris G Antonopoulos

Athanasios E Manos

Charalampos D Skokos

© IC-SCCE SALI: AN EFFICIENT INDICATOR OF CHAOS WITH APPLICATION TO 2 AND 3 DEGREES OF FREEDOM HAMILTONIAN SYSTEMS

Keywords: Hamiltonian systems, Chaotic motion, Ordered motion, Lyapunov Exponents

The Smaller Alignment Index (SALI) is a very useful and efficient indicator that can distinguish rapidly and with certainty between ordered and chaotic motion in Hamiltonian systems. This is based on the different behavior of the SALI in the two cases, as it fluctuates around non-zero values for ordered orbits and converges exponentially to zero for chaotic orbits. In this paper we present a detailed numerical study, comparing the advantages of the SALI to those of the Maximal Lyapunov Exponent (MLE), which is the traditional indicator of chaos detection in dynamical systems. Exploiting these advantages, we demonstrate how the SALI can be used to identify even tiny regions of order and chaos in Hamiltonian systems of 2 and 3 degrees of freedom. 1

INTRODUCTION

The distinction between ordered and chaotic motion in dynamical systems is fundamental in many areas of applied sciences. This distinction is particularly difficult in systems with many degrees of freedom (dof), basically because it is not feasible to visualize their phase space. So, we need fast and accurate tools to give us information about the chaotic or ordered character of orbits, especially for conservative systems.

Many methods have been developed over the years trying to give an answer to this problem. The inspection of the successive intersections of an orbit with a Poincaré surface of section (PSS) [START_REF] Lieberman | Regular and Chaotic Dynamics[END_REF] has been used mainly for 2 dimensional (2d) maps and 2 dof Hamiltonian systems. One of the most common methods of chaos detection is the computation of the maximal Lyapunov Characteristic Number (LCN) [START_REF] Benettin | Kolmogorov entropy and numerical experiments[END_REF][START_REF] Froeschlé | The Lyapunov characteristic exponents -Applications to celestial mechanics[END_REF] , which can be applied for systems with many degrees of freedom. Another very efficient method is the "Frequency Map Analysis" developed by Laskar [START_REF] Laskar | The chaotic motion of the solar system -A numerical estimate of the size of the chaotic zones[END_REF][START_REF] Laskar | The measure of chaos by the numerical analysis of the fundamental frequencies. Application to the standard mapping[END_REF] . In recent years new methods have been introduced such as the study of spectra of "Short Time Lyapunov Characteristic Numbers" [START_REF] Froeschlé | Generalized Lyapunov characteristic indicators and corresponding Kolmogorov like entropy of the standard mapping[END_REF][START_REF] Lohinger | Generalized Lyapunov exponents indicators in Hamiltonian dynamics -an application to a double star system[END_REF] or "Stretching Numbers" [START_REF] Voglis | Invariant spectra of orbits in dynamical systems[END_REF][START_REF] Contopoulos | Invariant spectra in Hamiltonian systems[END_REF] and the "Spectral Distance" of such spectra [START_REF] Voglis | Detection of ordered and chaotic motion using the dynamical spectra[END_REF] , as well as the study of spectra of helicity and twist angles [START_REF] Contopoulos | Spectra of stretching numbers and helicity angles in dynamical systems[END_REF][START_REF] Contopoulos | A fast method for distinguishing between ordered and chaotic orbits[END_REF][START_REF] Froeschlé | Twist angles: a method for distinguishing islands, tori and weak chaotic orbits. Comparison with other methods of analysis[END_REF] . In addition, Froeschlé C., et al. [START_REF] Froeschlé | Fast Lyapunov indicators. Application to asteroidal motion[END_REF] introduced the Fast Lyapunov Indicator (FLI), while Vozikis et al. [START_REF] Vozikis | The power spectrum of geodesic divergences as an early detector of chaotic motion[END_REF] proposed a method based on the frequency analysis of "Stretching Numbers".

Recently, a new, fast and easy to compute indicator of the chaotic or ordered nature of orbits, has been introduced [START_REF] Skokos | Alignment indices: a new, simple method for determining the ordered or chaotic nature of orbits[END_REF] and was applied successfully in different dynamical systems [START_REF] Skokos | Alignment indices: a new, simple method for determining the ordered or chaotic nature of orbits[END_REF][START_REF] Skokos | Smaller Alignment Index: Detecting order and chaos in conservative dynamical systems[END_REF][START_REF] Skokos | Smaller Alignment Index (SALI): Determining the ordered or chaotic nature of orbits in conservative dynamical systems[END_REF][START_REF] Skokos | How does Smaller Alignment Index (SALI) distinguish order from chaos?[END_REF][START_REF] Skokos | Detecting order and chaos in Hamiltonian systems by the SALI method[END_REF][START_REF] Panagopoulos | Existence and stability of localized oscillations in 1-dimensional lattices with soft spring and hard spring potentials[END_REF][START_REF] Voglis | Mass components in ordered and in chaotic motion in galactic N-body models[END_REF][START_REF] Széll | Investigation of the Caledonian Symmetric Four-Body Problem[END_REF] , the so-called "Smaller Alignment Index (SALI)". In the present paper, we first recall the definition of the Smaller Alignment Index and then we show its effectiveness in distinguishing between ordered and chaotic motion, by applying it to Hamiltonian systems of 2 and 3 dof.

DEFINITION OF THE SMALLER ALIGNMENT INDEX (SALI)

Let us consider the n -dimensional phase space of a conservative dynamical system, which could be a symplectic map or a Hamiltonian flow. We consider also an orbit in that space with initial condition 1 2 n P(0)=(x (0),x (0), ,x (0)) … and a deviation vector ξ(0)=(dx (0),dx (0), ,dx (0)) … from the initial point P(0) . In order to compute the SALI for a given orbit one has to follow the time evolution of the orbit with initial condition P(0) itself and also of two deviation vectors 1 2 ( ), ( ) t t ξ ξ which initially point in two different directions. The evolution of these deviation vectors is given by the variational equations for a flow and by the tangent map for a discrete-time system. At every time step the two deviation vectors 1 ( ) t ξ and 2 ( ) t ξ are normalized and the SALI is then computed as

1 2 1 2 1 2 1 2 ( ) ( ) ( ) ( ) ( ) min , ( ) ( ) ( ) ( ) t t t t SALI t t t t t ξ ξ ξ ξ ξ ξ ξ ξ ⎧ ⎫ ⎪ ⎪ = + - ⎨ ⎬ ⎪ ⎪ ⎩ ⎭ (1)
The properties of the time evolution of the SALI clearly distinguish between ordered and chaotic motion as follows: In the case of Hamiltonian flows or ndimensional symplectic maps with 2 n ≥ , the SALI fluctuates around a non-zero value for ordered orbits, while it tends to zero for chaotic orbits [START_REF] Skokos | Alignment indices: a new, simple method for determining the ordered or chaotic nature of orbits[END_REF][START_REF] Skokos | Smaller Alignment Index: Detecting order and chaos in conservative dynamical systems[END_REF][START_REF] Skokos | Smaller Alignment Index (SALI): Determining the ordered or chaotic nature of orbits in conservative dynamical systems[END_REF][START_REF] Skokos | How does Smaller Alignment Index (SALI) distinguish order from chaos?[END_REF][START_REF] Skokos | Detecting order and chaos in Hamiltonian systems by the SALI method[END_REF] . In the case of 2d maps the SALI tends to zero both for ordered and chaotic orbits, following however completely different time rates [START_REF] Skokos | Alignment indices: a new, simple method for determining the ordered or chaotic nature of orbits[END_REF] , which again allows us to distinguish between the two cases.

In 2 and 3 dof Hamiltonian systems the distinction between ordered and chaotic motion is easy because the ordered motion occurs on a 2d or 4d torus respectively on which any initial deviation vector becomes almost tangent after a short transient period. In general, two different initial deviation vectors become tangent to different directions on the torus, producing different sequences of vectors, so that SALI does not tend to zero but fluctuates around positive values [START_REF] Skokos | How does Smaller Alignment Index (SALI) distinguish order from chaos?[END_REF] . On the other hand, for chaotic orbits, any two initially different deviation vectors tend to coincide in the direction defined by the nearby unstable manifold and hence either coincides with each other, or become opposite. This means that the SALI tends to zero when the orbit is chaotic and to a nonzero value when the orbit is ordered. Thus, the completely different behavior of the SALI helps us to distinguish between ordered and chaotic motion in Hamiltonian systems with 2 and 3 dof and in general in dynamical systems of higher dimensionality.

APPLICATION OF THE SMALLER ALIGNMENT INDEX 3.1 The Hénon-Heiles 2 dof Hamiltonian system

In order to illustrate the effectiveness of the SALI method in the determination of the dynamical character of orbits (e.g. chaotic or ordered), we consider the 2 dof Hénon-Heiles Hamiltonian:

2 2 2 2 2 3 1 1 1 ( , , , ) ( ) ( ) 2 2 3 x y x y H x y p p p p x y x y y E = + + + + - ≡ (2)
where 1 8 E = is the constant total energy of the system, x , y are the generalized coordinates and x p , y p are the generalized conjugate momenta [START_REF] Hénon | The applicability of the third integral of motion: Some numerical experiments[END_REF] .

For the system (2) we consider three characteristic individual orbits and compute for each one the evolution of the SALI. The first orbit with initial condition ( , , , ) (0, 0.1, 0.490578,0)

x y

x y p p = we call it R, the second orbit with initial condition ( , , , ) (0, 0.0135, 0.499816,0)

x y

x y p p =we call it S and finally the orbit with initial condition (0, 0.015, 0.499773,0) we call it C. We also use 1 (0) (1, 0, 0, 0) ξ = and 2 (0) (0, 0,1, 0) ξ = as initial deviation vectors.

As we can see from figure 1(a) the SALI of the orbit R (solid line) fluctuates around SALI 1.25 ≈ remaining almost constant after [START_REF] Laskar | The measure of chaos by the numerical analysis of the fundamental frequencies. Application to the standard mapping[END_REF] 10 time steps indicating the ordered character of the orbit, while the SALI of orbit C (dash dot line) falls abruptly reaching the limit of the accuracy of the computer precision ( [START_REF] Skokos | Alignment indices: a new, simple method for determining the ordered or chaotic nature of orbits[END_REF] 10 -) after about 1000 time steps. This implies that the two deviation vectors are represented by the same or opposite sign numbers in the computer and thus they have been aligned with the direction of the most unstable nearby manifold. On the other hand, the SALI of orbit S (dot line) shows an intermediate behaviour. This orbit is called "sticky" as it remains at the borders of an island of stability for a long time but eventually enters the chaotic sea. Thus, the SALI shows the chaotic character of the orbit after about 5000 time steps when the "stickiness" phenomenon stops to play the dominant role in the dynamics. After this time, the SALI shows again the same behaviour as for the chaotic orbit C.

In figure 1(b) we see the evolution of the MLE of the same three orbits as a function of time (in logarithmic scale). The computation of MLE does not help us to decide with certainty the nature of the orbit as fast as the SALI. This becomes evident by comparing the evolution of the MLE and the SALI in figure 1. We can be absolutely certain that orbit S is chaotic at start to show a deviation from this decrease, although, in order to get a clear evidence of convergence to a non-zero value, further computations are needed. Thus, using the SALI method we are able to determine the nature of orbits with great certainty and in general earlier than the computation of the corresponding MLEs. In order to present the effectiveness of the SALI in detecting even tiny regions of chaos and order, we compute the SALI for a very large grid of equally distributed initial conditions on the PSS ( , ) y p y of Hamiltonian (2), assigning a different color to each initial condition according to the computed value of the SALI (we refer to this procedure as the full PSS scanning). In this way, we obtain a clear picture of the dynamics of the system on the PSS as we see in figure 2.

In particular, in figure 2(a) we compute the SALI of every individual initial condition on the grid for 2000 t = and assign an appropriate color according to the value of the index: If and the corresponding color not only to the initial condition but also to every subsequent intersection of the orbit with the PSS taking also into account the symmetry of the Hamiltonian (2) with respect to y p (we refer to this procedure as the fast PSS scanning). In this way, orbits with initial conditions on the grid to which a color has already been assigned, as they were intersection points with the PSS of previously computed orbits, are not computed again and so the whole process is speeded up significantly. In figure 2, we see that the differences of the two panels are almost not visible. We note that following the above procedure we are able to distinguish very tiny islands of stability located inside the chaotic sea.

Figure 2. Regions of different values of the SALI on the PSS 0 x = of the Hénon-Heiles system (2) for a grid of 500 500 × equally spaced initial conditions. In (a) the SALI was computed for every initial condition (full PSS scanning) while in (b) the final value of the SALI is attributed not only to the initial condition but also to all successive intersection points of the orbits with the PSS.

Chris G. Antonopoulos, Athanasios E. Manos and Charalampos D. Skokos By applying the above procedure and assuming that an orbit is chaotic beyond any doubt if 8 SALI 10 - < , we computed the percentage of the chaotic orbits in the whole PSS. Similar results can be found in the paper of Hénon and Heils [START_REF] Hénon | The applicability of the third integral of motion: Some numerical experiments[END_REF] where the authors used a different numerical technique in finding the percentage of chaotic orbits. The results of the full and fast PSS scanning are presented in figure 3. In figure 3(a), we see that both algorithms indicate that for energy 1 8 E = the PSS of figure 2 consists of about 55% of chaotic orbits as for 2000 t > both curves stabilize close to each other and further computations do not practically change the results. In figure 3(b) we plot the corresponding evolution of the CPU time on a Pentium 4 2GHz PC (in sec) needed for assigning a value of the SALI to every grid point on the PSS of system (2) as a function of the total time t needed to integrate all the individual orbits. Practically, this is the CPU time needed for obtaining the data of figure 2(b). The CPU time increases for small integration times t but after 2000 t exhibits a plateau as time t is long enough for assigning a SALI value to many grid points, that we do not take as new initial conditions, since the computed orbits have many intersections with the PSS. As t increases more, the CPU time slightly grows due to the increase of the integration time. So, we realize that although the fast PSS scanning is not absolutely accurate (the two curves of figure 3(a) do not coincide) it is preferred for the complete study of big regions of the PSS as the CPU time required is significantly less than the CPU time needed for the full PSS scanning. For example, the full PSS scanning algorithm needs about 68300 sec of CPU time to integrate all the initial conditions on the grid for 10000 t =

, while the fast PSS scanning algorithm needs only about 2100 sec!

The SALI method can be used also to discriminate between order and chaos, not only on the PSS of a system, but also in any parameter space. As an example, in figure 4, we plot the SALI values on the plane ( , ) y E for the Hamiltonian (2). In particular, we take orbits with initial conditions 0 x = , [ 0. . This is a similar plot to the one obtained by Fouchard et. al [START_REF] Fouchard | On the relationship between fast Lyapunov indicator and periodic orbits for continuous flows[END_REF] using the FLI method for chaos detection. From figure (4) we see that as energy E increases from 0.08 to 0.18 initial conditions with [ 0.3,0.4] y ∈that were initially ordered became chaotic, indicating that a lot of KAM tori (islands of stability) have been destroyed and in their position chaotic regions have been appeared. We stop our computations for 0.18 E = as for 0.18 E > escapes are possible.

The relation between the SALI and the Lyapunov exponents in a 3 dof Hamiltonian system

In order to illustrate the relation between the SALI and the 2 MLEs 1 L and 2 L respectively, we consider the 3 dof Hamiltonian system [START_REF] Contopoulos | On the number of isolating integrals in Hamiltonian systems[END_REF]

3 2 2 2 2 1 2 3 1 2 3 1 2 1 3 1 ( , , , , , ) ( ) 2 i i i i
H q q q p p p q p q q q q E ω where 1 2 , q q and 3 q are the generalized coordinates, 1 2 , p p and 3 p are the generalized conjugate momenta and 0.09 E = is the total constant energy of the system. Following [27], we use in our calculations for compatibility reasons. The characteristic property of system (3) according to [27] is that the system possesses two kinds of chaotic regions in its phase space. It has a large stochastic region where the two MLEs are of the same magnitude, called the "big sea" and stochastic regions, which are apparently separated from the "big sea" called "small seas". In these "small seas" the two MLEs are not of the same magnitude and they differ significantly.

= = + + + ≡ ∑ (3)
According to the recently published paper by Skokos et. al [START_REF] Skokos | Detecting order and chaos in Hamiltonian systems by the SALI method[END_REF] , the evolution of the SALI is related to the two largest Lyapunov exponents 1 L and 2 L of a 3 dof Hamiltonian system by (q ,q ,q ,p ,p ,p )=(0,0,0,0.14142,0.11892,0.29428) of the 3 dof Hamiltonian In order to support the validity of equation ( 4), we consider the orbit of system (3) with initial condition (q ,q ,q ,p ,p ,p )=(0,0,0,0.14142,0.11892,0.29428) originated at the "big stochastic sea". In figure 5(a) we plot the time evolution of the two largest Lyapunov exponents 1 ( ) L t and 2 ( ) L t of the orbit which are computed by

Figure 1 .

 1 Figure 1. Evolution of (a) the SALI and (b) the MLE with respect to time t in log-log scale for the orbits named R, S and C (see text).

  is colored light gray. In producing figure2(b), we performed a slightly different procedure to speed up the whole process without having significant errors: we assign the final value of the SALI at 2000 t =

Figure 3 .

 3 Figure 3. (a): The percentage of the chaotic orbits (having

Figure 4 :

 4 Figure 4: The SALI values for a grid of 500 500 × equally spaced initial conditions in the plane ( , ) y E for the Hénon-Heiles system (2) each integrated up to 2000 t =. If

Figure 5 .

 5 Figure 5. (a): The evolution of the two largest Lyapunov exponents 1 ( )L t and 2 ( ) L t for the chaotic orbit with initial condition 1 2 3 1 2 3 (q ,q ,q ,p ,p ,p )=(0,0,0,0.14142,0.11892,0.29428) of the 3 dof Hamiltonian (3). (b): Time

  Figure 5. (a): The evolution of the two largest Lyapunov exponents 1 ( )L t and 2 ( ) L t for the chaotic orbit with initial condition 1 2 3 1 2 3 (q ,q ,q ,p ,p ,p )=(0,0,0,0.14142,0.11892,0.29428) of the 3 dof Hamiltonian (3). (b): Time evolution of the SALI for the same orbit (black solid line) and the functions

ACKNOWLEDGMENTS Ch. Antonopoulos was partially supported by the "Karatheodory" graduate student fellowship No 2464 of the University of Patras and by the "Heraclitus" research program of the Greek Ministry of Development. Ch. Skokos was partially supported by the Research Committee of the Academy of Athens.

the method proposed by Benettin et. al. [START_REF] Benettin | Meccanica[END_REF] . We note that at final time 5 3 10 t ≈ ⋅ both exponents are of the same order of magnitude as 1 0.03747 L ≈ and 2 0.0105 L ≈ . In figure 5(b), we plot the SALI of the orbit as a function of time t (in linear scale) as well as the approximated functions and 3 c = . We see that using the SALI method we determine the chaotic nature of the orbit much faster than the computation of the Lyapunov exponents as for 600 t

we have 8 SALI 10 -while 1 ( ) L t and 2 ( ) L t still do not show clear evidence of convergence. We also note that the second Lyapunov exponent 2 L plays a crucial role in describing the behaviour of the SALI, as

CONCLUSIONS

In this paper we have illustrated the advantages of the SALI method in distinguishing between order and chaos in 2 and 3 dof Hamiltonian systems. Our results can be summarized as follows:

• The SALI proves to be an ideal indicator of chaoticity independent of the dimensions of the system. It tends to zero for chaotic orbits while it fluctuates around non-zero positive values for ordered ones, distinguishing clearly between the two cases.

• The SALI values characterize an orbit of being chaotic or ordered. Exploiting this feature of the index, we have plotted detailed phase space portraits for a 2 dof Hamiltonian system where the chaotic and ordered regions are clearly distinguished. We were thus able to trace very small islands of ordered motion, whose detection by traditional methods would be much more difficult and time consuming. • We have also presented an example of a chaotic orbit in a 3 dof Hamiltonian system, whose SALI tends to zero following an exponential rate which is related to the difference of the two largest Lyapunov exponents 1 L and 2 L . In this particular case, although 2 L is smaller than 1 L , both exponents have the same order of magnitude and so, 2

L cannot be neglected from the description of the time evolution of the SALI.


A Pilotto

P Dollfus

J Saint-Martin

M Pala

Full Quantum Simulation of Shockley-Read-Hall Recombination in p-i-n and Tunnel Diodes

Keywords: NEGF, SRH Recombination, Semiconductor Diodes

We have included Shockley-Read-Hall (SRH) generation/recombination in Non-equilibrium Green's function (NEGF) calculations via a multiphonon relaxation model. The model has been used to study how the presence of defects affects the current-voltage characteristics of GaAs p-i-n diodes, and of an InGaAs tunnel diode. Regarding p-i-n diodes, we show that SRH generation/recombination is responsible for ideality factors approaching the theoretical value of two in the forward bias regime, while in reverse bias the recombination current density varies slowly with the applied voltage. In all the considered cases, the defects located in the center of the active region proved to be the most effective in allowing trap-assisted tunneling from the valence to the conduction band. Finally, the inclusion of SRH recombination in NEGF simulations of Esaki tunnel diodes permits to predict a realistic degradation of the peak-to-valley current ratio.

Introduction

Trap assisted generation/recombination of electron and holes occurring in the active region of electron devices is known to be a limiting factor of their electrical performance. For instance, the dark count rate of single photon avalanche diodes (SPADs) is largely influenced by impact ionization of thermally generated carriers that, at large electric fields, trigger undesired avalanche multiplication and induce false photon detection [START_REF] Bronzi | Spad figures of merit for photon-counting, photon-timing, and imaging applications: A review[END_REF][START_REF] Panglosse | Dark count rate modeling in singlephoton avalanche diodes[END_REF]. At the same time, empty traps can be filled during the avalanche subsequent to the generation of electron-hole pairs induced by photon absorption. Detrapping of these charged states at a later time is the cause of secondary pulses in the current waveform (after pulsing) that are erroneously interpreted as photon counts [START_REF] Bronzi | Spad figures of merit for photon-counting, photon-timing, and imaging applications: A review[END_REF][START_REF] Vilà | Geiger-mode avalanche photodiodes in standard CMOS technologies[END_REF]. On the other hand, in Esaki diodes, recombination induced by defect states is one of the main contributors to the degradation of the peak-to-valley current ratio in the forward bias regime [START_REF] Meyerhofer | Degenerate Germanium. I. tunnel, excess, and thermal current in tunnel diodes[END_REF][START_REF] Schenk | Tunneling between density-of-state tails: Theory and effect on Esaki diodes[END_REF].

An accurate modeling of defect-mediated generation/recombination is therefore of fundamental importance for the correct assessment of the behavior of semiconductor diodes and, in the literature, many models to treat this problem can be found, ranging from semiclassical [START_REF]Sentaurus Device User Guide[END_REF] to full-quantum approaches [START_REF] Pala | Interface traps in InAs nanowire tunnel-FETs and MOSFETs-part I: Model description and single trap analysis in tunnel-FETs[END_REF][START_REF] Shao | Sub-10-nm diameter GaSb/InAs vertical nanowire Esaki diodes with ideal scaling behavior: Experiments and simulations[END_REF][START_REF] Shao | Scaling of GaSb/InAs vertical nanowire Esaki diodes down to sub-10-nm diameter[END_REF]. Recently, a model to include Shockey-Read-Hall (SRH) type recombination [START_REF] Shockley | Statistics of the recombinations of holes and electrons[END_REF][START_REF] Hall | Electron-hole recombination in Germanium[END_REF] via nonradiative multiphonon relaxation [START_REF] Schenk | An improved approach to the Shockley-Read-Hall recombination in inhomogeneous fields of space-charge regions[END_REF] in NEGF calculations has been proposed in [START_REF] Aeberhard | Nonequilibrium Green's function picture of nonradiative recombination of the Shockley-Read-Hall type[END_REF], and successfully applied to study forward biased p-i-n junctions [START_REF] Aeberhard | Nonequilibrium Green's function picture of nonradiative recombination of the Shockley-Read-Hall type[END_REF] and LEDs [START_REF] Gonzalez Montoya | Nonequilibrium Green's function modeling of trap-assisted tunneling in In x Ga 1-x N/GaN Light-Emitting Diodes[END_REF]. In this paper, we use the model of [START_REF] Aeberhard | Nonequilibrium Green's function picture of nonradiative recombination of the Shockley-Read-Hall type[END_REF] to study the impact of SRH recombination occurring at different device regions on the forward and reverse current-voltage characteristics (I -V ) of GaAs p-i-n junctions and of an InGaAs Esaki diode. The paper proceeds as follows: Section 2 briefly describes the model of [START_REF] Aeberhard | Nonequilibrium Green's function picture of nonradiative recombination of the Shockley-Read-Hall type[END_REF] and its implementation. Results obtained by simulations of bulk GaAs, GaAs p-i-n diodes, and of an InGaAs Esaki diode are shown in Section 3. Finally, Section 4 reports our concluding remarks.

Model Description

The model for the treatment of SRH recombination via multiphonon relaxation in NEGF calculations [START_REF] Aeberhard | Nonequilibrium Green's function picture of nonradiative recombination of the Shockley-Read-Hall type[END_REF] assumes that the defect density is strongly localized in space and energy, and, in the case of a structure with confinement along the z direction, it is therefore represented by a Dirac's delta function at z = z def and E = E def .

Here, as in [START_REF] Aeberhard | Nonequilibrium Green's function picture of nonradiative recombination of the Shockley-Read-Hall type[END_REF][START_REF] Gonzalez Montoya | Nonequilibrium Green's function modeling of trap-assisted tunneling in In x Ga 1-x N/GaN Light-Emitting Diodes[END_REF], the Dyson's equation for the retarded Green's function of the defect (G d ) is not solved, but, instead, the lesser-than and greaterthan Green's functions (G < d and G > d , respectively) are computed by using the quasi-equilibrium expressions

G < d (E) = if d (E)A d (E), (1)
G > d (E) = -i(1 -f d (E))A d (E), (2)
where A d is the defect's spectral function and f d describes the defect's occupation. The value of f d (E = E def ) can be derived from the balance at steady-state between the capture and emission rates [START_REF] Aeberhard | Nonequilibrium Green's function picture of nonradiative recombination of the Shockley-Read-Hall type[END_REF], namely

f d (E def ) = -i Σ < dc (E def ) + Σ < dv (E def ) Γ dc (E def ) + Γ dv (E def ) , (3)
where Γ = i(Σ > -Σ < ).

In the case of electrons, the self-energies that describe carrier capture/emission into/from the defect state and into/from the bands are computed with Eqs. ( 4) and ( 5), respectively.

Σ ≶ dc (z def , E) = 1 A k l≥0 M em/capt (l) dzG ≶ c (k , z, z, E + l Ω 0 )δ(z -z def ) (4)
Σ ≶ cd (z, z, E) = ρ d l≥0 M em/capt (l)G ≶ d (E -l Ω 0 )δ(z -z def ) (5)
M capt (l) = M 0 dc (l -S) 2 S e -S(2f BE +1) f BE + 1 f BE l/2 I l (2S f BE (f BE + 1)) (6)
where ρ d is the defect density, Ω 0 is the phonon frequency, l identifies the number of phonons involved in the recombination process, S is the Huang-Rhys factor, f BE is the phonon occupation given by the Bose-Einstein distribution, I l is the modified Bessel function of l-th order, and M em (l) =

M capt e -l Ω 0 /k B T . Expressions similar to Eqs. ( 4), ( 5) are used to compute the self-energies (Σ dc and Σ vd ) for the coupling between the defect and valence band states. The constants M 0 dc and M 0 dv are linked to the carrier thermal velocities (v th,n/p ) and capture cross sections (σ n/p ) as described in [START_REF] Aeberhard | Nonequilibrium Green's function picture of nonradiative recombination of the Shockley-Read-Hall type[END_REF].

Finally, current conservation is obtained by self-consistently solving the equations for the retarded and the lesser(greater)-than Green's functions.

Results

This section reports the results of the simulations of bulk GaAs, GaAs p-i-n diodes and an InGaAs p-i-n diode when the method reported in Sec. 2 is used. In the following, the parameters reported in Tab. 1 have been employed.

Defect states have always assumed to be located at midgap.

Material m e [m 0 ] m h [m 0 ] E g [eV] S Ω 0 [

Bulk GaAs

As a validation of our implementation, we have verified that the recombination rate computed by using the aforementioned approach (Eq. ( 7)) is equivalent to the semiclassical expression (Eq. ( 8)) [START_REF] Shockley | Statistics of the recombinations of holes and electrons[END_REF] when the the carrier lifetimes are τ n/p = ρ d σ n/p v th,n/p and the electron and hole densities are obtained from NEGF calculations. The results in Fig. 1, obtained by using an effective mass Hamiltonian with two parabolic bands (see Tab. 1), and for a defect density ρ d = 10 14 cm -3 show an excellent agreement between Eqs. [START_REF] Pala | Interface traps in InAs nanowire tunnel-FETs and MOSFETs-part I: Model description and single trap analysis in tunnel-FETs[END_REF] and [START_REF] Shao | Sub-10-nm diameter GaSb/InAs vertical nanowire Esaki diodes with ideal scaling behavior: Experiments and simulations[END_REF]. Results obtained by using the NEGF formalism (Eq. ( 7), symbols) are compared to the semiclassical expression [START_REF] Shockley | Statistics of the recombinations of holes and electrons[END_REF] (Eq. ( 8), lines) for σ n =σ p =10 -14 cm 2 and ρ d =10 14 cm -3 .

U SRH = ρ d dE 2π [Σ < dc (E)G > d (E) -Σ > dc (E)G < d (E)] (7)
U SRH = np -n eq p eq τ p (n + n 1 ) + τ n (p + p 1 ) (8

GaAs p-i-n Diodes

We factors that approach or exceed two (see Fig. 3). Ideality factors larger than two arise from the enhancement of the trap-assisted tunneling process pro-vided by the band bending; this phenomenon is particularly relevant for thin diodes. On the other hand, in reverse bias configurations the recombination current is slowly varying with the applied bias, except for those cases where a single defect is located far from the center of the intrinsic region. In fact, by looking at the spectral current densities in Fig. 4, we notice that, due to the position of the defect with respect to the Fermi levels of the contacts, tunneling mediated by traps located at the center of the intrinsic region is the one that mostly affects the recombination current, while recombination at other locations becomes significant only for large applied bias voltages.

J rec [A m -2 eV -1 ] J rec [A m -2 eV -1 ] μ p μ p μ p μ p μ n μ n μ n μ n a) b) c) d) V = 0.48 V V = 0.84 V V = -0.6 V V = -1.2 V

Conclusions

We

) 2 Figure 1 :

 21 Figure 1: Recombination rate in bulk GaAs as a function of the applied bias voltage.

  have employed the model to simulate the recombination current characteristics of GaAs p-i-n diodes with different thickness of the intrinsic region (L = 40 nm and L = 80 nm). A doping concentration N d = N a = 10 17 cm -3 has been assumed in both contact regions and a linear potential profile across the intrinsic layer has been considered. For the two diodes, five different traps configurations have been analyzed: a single midgap defect located at z def = L/4, L/2, or 3L/4 and two random distributions of five traps, namely R1 40 = [4 nm, 7.6 nm, 10 nm, 23.3 nm, 39.2 nm], R2 40 = [10 nm, 15 nm, 20 nm, 25 nm, 30 nm] for the 40 nm diode, and R1 80 = [10 nm, 11.6 nm, 15.2 nm, 32.8 nm, 57.2 nm], R2 80 = [20 nm, 30 nm, 40 nm, 50 nm, 60 nm] for the thicker diode. Recombination currents are shown in Fig. 2. As expected, in forward bias the recombination current increases with ideality

Figure 2 :

 2 Figure 2: Recombination current as a function of the applied voltage for GaAs p-i-n diodes with intrinsic regions thickness (a) L=40 nm and (b) 80 nm. Different configurations have been considered: single midgap defects located at L/4 (black circles), L/2 (red triangles), 3L/4 (green squares) and two random distributions of five traps (blue plus signs and magenta crosses). σ n =σ p =5 × 10 -14 cm 2 and ρ d =10 13 cm -3 .

Figure 3 :

 3 Figure 3: Ideality factors in forward bias for the p-i-n diodes and the trap configurations of Fig. 2.

Figure 4 :Figure 5 :Figure 6 :

 456 Figure 4: Spectral recombination current density at (a) V =0.48 V, (b) V =0.84 V, (c) V = -0.6 V, and (d) V = -1.2 V of a 40 nm long GaAs p-i-n diode when five midgap defects from z=10 nm to z=30 nm are simulated (blue circles). σ n =σ p =5 × 10 -14 cm 2 and ρ d =10 13 cm -3 .

  have employed a quantum model to include SRH recombination in NEGF simulations. The model has been used to study how the defect location in the active region of p-i-n diodes affects the characteristics of these devices, showing that SRH generation/recombination via multiphonon relaxation mostly occurs at defects located close to the center of the diode's active region, due to the relative position of these trap states with respect to the Fermi levels of the contacts. In the simulation of Esaki tunnel diodes, the inclusion of the SRH recombination is able to predict the increase of the valley currents and the corresponding degradation of the peak-to-valley current ratio (from 3400 to 22). Therefore, this model can shed light on the role played by traps in the current characteristics of low-power devices such as the tunnel-FET.

Table 1 :

 1 Parameters used throughout this paper.

	meV]

Acknowledgments

This work was supported by the ANR project "GeSPAD" (ANR-20-CE24-0004).


Mélanie Miranda

Daan Noordermeer

Benoit Moindrot
email: benoit.moindrot@i2bc.paris-saclay.fr

Detection of allele-specific 3D chromatin interactions using high-resolution innucleus 4C-seq

Keywords: Circular Chromosome Conformation Capture (4C), Nuclear organization, DNA loops, Topologically associating domains (TADs), High-throughput sequencing, Singlenucleotide polymorphism (SNP)

Chromosome Conformation Capture techniques are a set of methods used to determine 3D genome organization through the capture and identification of physical contacts between pairs of genomic loci. Among them, 4C-seq (Circular Chromosome Conformation Capture coupled to high-throughput sequencing) allows the identification and quantification of the sequences interacting with a pre-selected locus of interest. 4C-seq has been widely used in the literature, mainly to study chromatin loops between enhancers and promoters or between CTCF binding sites and to identify chromatin domain boundaries.

As 3D-contacts may be established in an allele-specific manner, we describe an upto-date allele-specific 4C-seq protocol, starting from the selection of allele-specific viewpoints to Illumina sequencing. This protocol has mainly been optimized for cultured mammalian cells, but can be adapted for other cell types with relatively minor changes in initial steps.

Introduction

The 3D organization of genomes plays important roles in the transcriptional regulation of genes [1]. Over the past two decades, the development of "chromosome conformation capture" assays has deepened our understanding of how chromosomes are folded, and the regulatory consequences of this folding [2,3]. In particular, the discovery of the relatively insulated Topologically Associating Domains (TADs) shed light on an important regulatory level [START_REF] Nora | Spatial partitioning of the regulatory landscape of the Xinactivation centre[END_REF][START_REF] Dixon | Topological domains in mammalian genomes identified by analysis of chromatin interactions[END_REF][6][7][8]. The formation of the TADs is largely driven by the coordinated action of CTCF and the Cohesin complex [9][START_REF] Fudenberg | Formation of Chromosomal Domains by Loop Extrusion[END_REF][START_REF] Chang | TADs and Their Borders: Free Movement or Building a Wall?[END_REF]. These two molecular actors are found enriched at the TAD boundaries, thereby restricting, according to current models, promoterenhancer communication to a single TAD [START_REF] Fudenberg | Formation of Chromosomal Domains by Loop Extrusion[END_REF][START_REF] Razin | Gene functioning and storage within a folded genome[END_REF][START_REF] Rao | A 3D map of the human genome at kilobase resolution reveals principles of chromatin looping[END_REF].

High-throughput sequencing-based Chromosome Conformation Capture (3C) assays come in different flavors [3]. Hi-C scores for interactions between all possible pairs of genomic loci at a resolution that is largely dictated by the sequencing depth [START_REF] Lieberman-Aiden | Comprehensive mapping of long-range interactions reveals folding principles of the human genome[END_REF]. 4Cseq (Circular Chromosome Conformation Capture) on the other hand queries for the chromatin interactions between a locus of interest (thereafter named the "viewpoint") and essentially any other genomic loci [START_REF] Noordermeer | The dynamic architecture of Hox gene clusters[END_REF][START_REF] Werken Hjg Van De, Landan | Robust 4C-seq data analysis to screen for regulatory DNA interactions[END_REF]. The resolution of 4C-seq assays is not limited by the sequencing power, but rather determined by the first restriction enzyme used. Capture-C presents an alternative way to determine the interaction frequencies between a relatively large set of genomic loci and the rest of the genome [START_REF] Hughes | Analysis of hundreds of cisregulatory landscapes at high resolution in a single, high-throughput experiment[END_REF][START_REF] Davies | Multiplexed analysis of chromosome conformation at vastly improved sensitivity[END_REF]. Yet, 4C-seq assays will be more flexible when interested in a few viewpoints (up to around 10 viewpoints), especially if 3D-contacts need to be described in an allele-specific manner.

The central principle of 3C is the detection of 3D interactions through the identification and the quantification of a ligation junction between two restriction fragments as a measure for initial proximity. High-resolution 4C-seq experiments start by the fixation of chromatin contacts using chemical crosslinking, for which formaldehyde is most commonly used. The chromatin is then digested using a restriction enzyme that recognizes a 4-bp sequence, followed by proximity ligation to create concatemers of restriction fragments whose linear arrangement directly reflects their initial spatial proximity. To characterize the ligation junctions between fragments, the resulting DNA is purified and subjected to a second round of enzymatic digestion and ligation, which results in circularization of the DNA. The circles containing the region of interest (viewpoint) and potential captured partners are amplified using an inverse PCR that concurrently generates ready-to-sequence libraries (Figure 1). The 3D partners are finally identified and quantified after Illumina sequencing of the ligation junctions and bioinformatical analysis of the resulting reads.

Although for most genomic intervals, the two alleles of a diploid cell are assumed to adopt a virtually identical average 3D conformation, there are a variety of cases where this may not be the case, for example in the context of X-chromosome inactivation [START_REF] Splinter | The inactive X chromosome adopts a unique three-dimensional conformation that is dependent on Xist RNA[END_REF], genomic imprinting [START_REF] Llères | CTCF modulates allele-specific sub-TAD organization and imprinted gene activity at the mouse Dlk1-Dio3 and Igf2-H19 domains[END_REF][START_REF] Noordermeer | Differential 3D chromatin organization and gene activity in genomic imprinting[END_REF] or cases of allelic exclusion [START_REF] Holwerda | Allelic exclusion of the immunoglobulin heavy chain locus is independent of its nuclear localization in mature B cells[END_REF][START_REF] Johanson | Genome organization in immune cells: unique challenges[END_REF][START_REF] Clowney | Nuclear aggregation of olfactory receptor genes governs their monogenic expression[END_REF]. Of particular relevance of this chapter, different strategies have been adopted to set-up allelic 4Cseq assays and determine allele-specific chromosomal contacts (Figure 2). These strategies fall into two distinct approaches. A first approach relies on an allele-specific restriction site, thereby solely amplifying the partners captured by one of the two alleles (Figure 2A and 2B). A second approach relies on the sequencing of a single-nucleotide variant located between one of the primers and a digested restriction site (Figure 2C and 2D). Both approaches have their own set of advantages and limitations (Figure 2). For instance, in the first approach, the design of viewpoints is less flexible but there is a much-reduced chance that reads are incorrectly attributed to the wrong allele, as restriction enzymes are very efficiently blocked by changes to the genomic sequence.

In contrast, the second approach is much easier for viewpoint design but there is a considerable chance that reads are attributed to the incorrect allele because the allelic identification is based on the sequencing of usually a single polymorphism.

In this chapter, we provide a detailed protocol of high-resolution and allele-specific 4Cseq assay that relies on an allele-specific restriction site. This is a further development from our previous protocol [START_REF] Matelot | Determination of High-Resolution 3D Chromatin Organization Using Circular Chromosome Conformation Capture (4C-seq)[END_REF], which includes a set of further improvements and adaptations. In particular, this in-nucleus version of the 4C-seq protocol incorporates the technical optimizations introduced for Hi-C experiments to preserve the nucleus prior ligation and reduce the frequency of random ligations occurring in diluted conditions [START_REF] Rao | A 3D map of the human genome at kilobase resolution reveals principles of chromatin looping[END_REF][START_REF] Nagano | Comparison of Hi-C results using in-solution versus in-nucleus ligation[END_REF]. Although a special focus is given to the design of allele-specific experiments, the following methods and guidelines will be useful for setting-up both allelic and non-allelic 4C-seq assays.

Materials

Equipment

1. Cooled centrifuges for 1.5 ml micro tubes and 15 ml and 50 ml conical tubes 2. Tube roller or mixer 3. Thermomixer with heating and cooling capability and holders for 1.5 ml micro tubes and 50 ml conical tubes. Alternatively, use a water bath.

4. Spectrophotometer capable of determining DNA concentration using dyeincorporation (e.g., QuBit Fluorometer, Thermo Fisher Scientific) 5. Magnetic separation rack for 1.5 ml micro tubes 6. PCR machine 7. -80°C freezer / -20°C freezer / Fridge 8. 15 ml and 50 ml conical tubes, 1.5 ml Safe-Lock tubes (e.g. Eppendorf), DNA 1.5 ml LoBind tubes (e.g. Eppendorf), 0.2 ml thin walled PCR tubes. 9. Tubes or cuvettes for DNA dye-incorporation (e.g. QuBit assay tubes)

Kits

1. Qiagen PCR Purification Kit (Qiagen, # 28104). If needed, additional PB buffer can be ordered separately (Qiagen, #19066)

Design of allele-specific 4C-seq primers

The construction of 4C-seq primers is indicated in Figure 3. From the 5'-end to the 3'end, they are composed of an Illumina P5 or P7 sequence for cluster generation, an Illumina Read Sequencing Primer sequence, and a viewpoint specific sequence that anneals to the 4C template. A variable index for multiplexing strategies can be added to the primer facing the restriction site of the second cutter (Figure 3).

For the allele-specific 4C-seq described here, the P5-containing primer faces the allele-specific restriction site of the first cutter (Figure 2B). This requires to know beforehand the location of allelic variants. In this respect, the mouse model represents a convenient system as a comprehensive list of SNPs in numerous mouse strains is available on the Mouse Genomes Project (see Note 2). This list needs to be intersected with the location of restriction sites to list allele-specific restriction sites.

The primer facing the allelic cut site of the first restriction enzyme should anneal as close as possible to restriction site (ideally within 30 bp), as Illumina reads will start from the 5'-end of the viewpoint specific sequence of the 4C-seq primer.

Each 4C primer pair should be tested and validated beforehand (see Note 3). 4C-seq primers with Illumina adapters are long and therefore needs to be PAGE-purified to eliminate shorter by-products (see Note 4). Primers are diluted to 100 μM and can be stored at -20°C.

Methods

See Notes 1 to 5 for important guidance regarding the experimental design before starting the experiment.

Cell crosslinking

1. Collect 10-15 x 10 6 cells and centrifuge (5 minutes, 300 g). Resuspend the cell pellet in 10 ml of Cell isolation buffer (see Note 1 if working with non-mammalian cells).

2. Centrifuge (5 minutes, 300 g, 4°C) and discard supernatant.

Resuspend the pellet in 12 ml of Crosslink buffer (2% formaldehyde final).

Gently pipet up-and-down to dissociate any aggregates. Incubate at room temperature for exactly 10 minutes on a tube roller mixer. Then, immediately transfer the tube on ice. 8. Centrifuge (1 minute, 300 g, 4°C) and discard supernatant. Snap-freeze the cell pellet in liquid nitrogen. Store the cell pellets at -80°C for up to 2 years.

Cell lysis and nucleus permeabilization

1. Gently resuspend the cell pellet in 1 ml of cold Cell lysis buffer. Dissociate any aggregates by pipetting up and down several times.

2. Transfer cells in 15 ml cold conical tube. Add 9 ml of cold Cell lysis buffer and incubate 20 minutes on ice.

3. Centrifuge the lysed cells (5 minutes, 2500 g, 4°C) and discard the supernatant.

4. Resuspend the pellet in 10 ml of cold Cell lysis buffer and incubate again 10 minutes on ice.

5. Centrifuge the lysed cells (5 minutes, 2500 g, 4°C) and discard the supernatant.

6. Resuspend the pellet with 500 μl of cold Cell lysis buffer and transfer the suspension to a 1.5 ml Safe-Lock tube.

7. Centrifuge (5 minutes, 2500 g, 4°C) and discard the supernatant. 9. Resuspend nuclei in 100 μl 0.5% SDS solution and incubate for exactly 10 minutes at 62°C while shaking in a Thermomixer (900 rpm).

First digestion

1. Add x µL of water and 25 µl of 20% Triton-100X. The volume x is calculated as follows: 323 µl minus the volume of enzyme used at step 3.3.5. This allows a final volume of 500 µl for the digestion.

2. Incubate for 15 minutes at 37°C while shaking in a Thermomixer (900 rpm).

See Note 6 if aggregates form at this or later steps.

3. Add 50 μl of the appropriate 10X enzyme restriction buffer and 2 μl of 20 mg/ml BSA.

4. Optional: Save a 5 μl aliquot of undigested material. Keep at 4°C or at -20°C.

5. Add 400 Units of the selected first restriction enzyme. The volume is now 500 µl. Mix well and incubate for 4 h at 37°C while shaking in a Thermomixer (900 rpm).

6. Add another 400 Units of the selected first restriction enzyme. Mix well and incubate overnight at 37°C while shaking in a Thermomixer (900 rpm).

DNA digestion efficiency

1. Save a 5-10 μl aliquot of digested chromatin in a 1.5 ml micro tube. Store the rest of the chromatin material at 4°C. Thaw the optional undigested material control (step 3.3.4).

2. Add 95-90 µl of water and 5 µl 5 M NaCl to both samples (final volume: 100 µL). 5. Add 2 μl of 10 mg/ml RNase A and incubate for 30 minutes at 37°C while shaking in a Thermomixer (750 rpm).

1. Add 120 μl phenol-chloroform-IAA to the solution. Shake vigorously and centrifuge (15 minutes, 16,000 g, room temperature). Phenol is toxic: work in a chemical fume hood with appropriate protection. 8. Wash the precipitate with a 70 % ethanol and centrifuge again (10 minutes, 16,000 g, 4°C) and carefully discard the supernatant. 9. Air-dry the precipitate and resuspend it in 20 μl of 10 mM Tris pH 8. 10. Visualize DNA on 1.5% agarose gel to evaluate digestion efficiency. Well digested material runs as a smear with DNA fragments ranging from 100 bp to a few kb. A light band of poorly-digested material can remain, but its size is usually slight smaller than undigested control (Figure 4A). If sample is underdigested (Figure 4A), add 400 U of the selected first restriction enzyme to the chromatin material, and incubate overnight at 37°C (900 rpm). If after 2 rounds of overnight digestion by the first restriction enzyme the chromatin remains under-digested, it is usually associated with poor cell lysis/permeabilization at initial steps. In this case, the cell lysis should be optimized on a new sample.

Restriction enzyme inactivation and First ligation

1. Take the chromatin sample out of the fridge and inactivate the first restriction enzyme by incubating the chromatin material at 65°C for exactly 20 minutes while shaking in a Thermomixer (600 rpm).

2. Prepare the 1.5 X ligation master mix as follows: 150 µl commercial 10X Ligation Buffer, 757 µl H2O, 75 µL 20% Triton-X100, 8 µl 20 mg/ml BSA.

3. Add up to 990 μl of ligation master mix for a final volume of 1.5 ml. Incubate at 16°C for 15 minutes while shaking in a Thermomixer (900 rpm).

4. Add 5 μl of T4 DNA ligase (20 U/μl, Weiss units). Incubate for 4h at 16°C while shaking in a Thermomixer (900 rpm).

De-crosslink samples by incubating overnight at 65°C while shaking in a

Thermomixer (900 rpm).

DNA Purification

2. Divide the samples over two (~750 μl) new Safe-Lock tubes and process both tubes as indicated in following steps. 12. Air-dry the precipitate for up to 1 hour at 37-55°C. [START_REF] Rao | A 3D map of the human genome at kilobase resolution reveals principles of chromatin looping[END_REF]. Dissolve the precipitate in 200 μl of ultrapure water. Incubate 1 hour at 37°C. This constitutes the 3C material.

14. Quantify a 2 µL aliquot of the 3C material using a DNA dye-incorporation assay for high concentration DNA (e.g., QuBit dsDNA BR assay kit; see Note 7).

Visualize 200 ng of 3C material on a 1% agarose gel to verify the ligation efficiency. Correctly ligated material will run at a considerably increased average length compared to digested control (highest average signal from several kb's up to over 10 kb).

Second digestion

1. Dilute the 3C material DNA to 100 ng/μl in the appropriate 1X restriction buffer for the second restriction enzyme (see Note 1 for the choice of restriction enzymes).

2. Add 1 Unit of the second restriction enzyme per μg of 3C material and incubate overnight at 37 °C (750 rpm).

Restriction enzyme inactivation

1. Inactivate the second restriction enzyme by incubating for 20 minutes at 65°C in a Thermomixer (600 rpm).

2. Take a 5 μl aliquot as 'digested control' and visualize on a 2% agarose gel.

Digested material should run as a smear ranging from 100 bp to 2 kb. If digestion is not complete, do another round of overnight digestion adding the same amount of the second restriction enzyme, and check the digestion efficiency the following day before continuing. 12. Centrifuge (1 minute, 16,000 g, room temperature) and discard the flowthrough.

13. Wash the column with 600 μl of Qiagen PE buffer. Centrifuge (1 minute, 16,000 g, room temperature).

14. Wash the column again with 600 µL of Qiagen PE buffer. Centrifuge (1 minute, 16,000 g, room temperature).

15. Transfer the 4 columns to 4 new 1.5 ml Safe-Lock tubes. [START_REF] Werken Hjg Van De, Landan | Robust 4C-seq data analysis to screen for regulatory DNA interactions[END_REF]. Eluted the 4C material by adding 50 μl of 10 mM Tris-HCl pH 7.5. Centrifuge (1 minute, 16,000 g, room temperature). 17. Combine the 4 elutions into a single 1.5 ml micro tube (total volume is around 180 µL).

18. Quantify a 2 µL aliquot of the 3C material using a DNA dye-incorporation assay for high concentration DNA (e.g. QuBit dsDNA BR assay kit; see Note 7). When starting from 10 7 mammalian cells, one should expect around 30 µg of 4C material at this step.

19. Store the 4C material at -20°C.

4C-seq Test PCR

Testing the correct amplification of 4C material should be conducted by a titration PCR using, ideally, a validated pair of 4C primers (see Note 3). The setup of the four PCR reactions is as follows:

Preparation of 4C-seq Illumina sequencing material

1. Do as many PCR reaction as necessary to amplify 1 -1.2 µg 4C material, using the highest concentration of 4C material that still provides linear amplification in the 4C-seq Test PCR. We routinely use 75 ng per PCR reaction and run a total of 16 PCR reactions (see Note 4 for information on Illumina primers and Note 8 for the importance of running multiple PCR reactions to generate sufficiently complex 4C material).

2. Pool the PCR reactions. Mix well and transfer a 100 µL aliquot into a new 1.5 ml tube. Amplicons in this aliquot will be purified SPRI beads in the following steps (see Note 9 for an alternative purification strategy). Store the remainder of the unpurified pooled PCR reactions at -20°C.

3. Allow the SPRI beads to warm-up to room temperature (30 min).

4. Vortex thoroughly the SPRI beads and immediately after add precisely 80 µL of the beads to the 100 µL aliquot from the PCR amplification (0.8 volumes beads to 1 volume of PCR reaction). Pipet up-and-down at least 10 times or vortex. Do a short spin to collect all liquid in the solution and incubate 5 minutes at room temperature.

5. Place the 1.5 ml microcentrifuge tube on a magnetic rack. Remove the supernatant once the solution is clear.

6. Gently add 1 ml freshly prepared 80% Ethanol while the tube remains on magnet, without disturbing the beads. Remove the Ethanol after 30 seconds.

7. Gently add again 1 ml 80% Ethanol while the tube remains on magnet, without disturbing the beads. Remove the Ethanol after 30 seconds. 8. Briefly spin the tube to pellet the remaining drops of Ethanol. Place the tube on the magnetic rack and carefully remove the residual ethanol. 9. Air-dry for 30 seconds. 12. Quantify a 4 µL aliquot of the PCR material using a DNA dye-incorporation assay for high concentration DNA (e.g., QuBit dsDNA BR assay kit).

13. Run 15 µL or up to 400 ng of purified 4C-seq PCR on a 1.5% Agarose gel.

Verify the correct amplification (see Note 3) and check for the absence of primer dimers.

14. The amplified material is ready to be sequenced. Multiplex by combining equal amount of 4C-seq Illumina sequencing material in a DNA LoBind Tube (see Note 10). A 3.3 ng/µL DNA concentration is equivalent to 20 nM when DNA fragments are on average 250 bp. Confirm the final concentration of the pooled 4C PCR material by quantifying a 2 µL aliquot using a DNA dye-incorporation assay for low concentration DNA (e.g. QuBit DNA HS assay kit). This material can be used directly for sequencing on an Illumina NextSeq 500 platform in a 75-bp Single-End run (see Note 11).

Data analysis

After Illumina sequencing, the resulting data should be analyzed using bioinformatics tools to create 4C-seq interaction profiles that depict the interaction frequency of each viewpoint with the rest of genome. Different 4C-seq analysis pipelines are available in the literature, which include 4Cseqpipe [START_REF] Werken Hjg Van De, Landan | Robust 4C-seq data analysis to screen for regulatory DNA interactions[END_REF], r3Cseq [START_REF] Thongjuea | r3Cseq: an R/Bioconductor package for the discovery of long-range genomic interactions from chromosome conformation capture and next-generation sequencing data[END_REF], FourCseq [START_REF] Klein | FourCSeq: analysis of 4C sequencing data[END_REF], Basic4Cseq [START_REF] Walter | Basic4Cseq: an R/Bioconductor package for analyzing 4C-seq data[END_REF], fourSig [START_REF] Williams | fourSig: a method for determining chromosomal interactions in 4C-Seq data[END_REF], w4Cseq [START_REF] Cai | w4CSeq: software and web application to analyze 4C-seq data[END_REF] or pipe4C [START_REF] Krijger | 4C-seq from beginning to end: A detailed protocol for sample preparation and data analysis[END_REF]. Another option, which we frequently used, was found in the recently defunct HTSstation [START_REF] David | HTSstation: a web application and open-access libraries for high-throughput sequencing data analysis[END_REF]. A stand-alone version of this analysis tool for 4C-seq data, combining python, R and bash scripts, is available upon request.

All 4C-seq data analysis pipelines roughly follow a similar data analysis strategy, including mapping of data, translation to interaction frequencies and either the identification of genomic regions with significantly enriched signal or the possibility to quantitatively compare between samples. Our preferred data analysis approach can be subdivided into five essential steps:

1. Demultiplexing and trimming of the reads. This is done by matching the 5'-end of the reads with the sequence of the 4C primers used for amplification, with some mismatches tolerated. Reads are trimmed at the location of the restriction sites for the following mapping step, with the part complementary to the primer sequence being removed.

Mapping of the reads to the reference genome using an alignment tool like

Bowtie2 [START_REF] Langmead | Fast gapped-read alignment with Bowtie 2[END_REF].

Virtual digestion of the reference genome to generate a list of valid restriction

fragments. This list is filtered out for restriction fragments lacking a cut site for the second restriction enzyme or restriction fragments containing repeated sequences at both extremities.

4. Assignment of read density to restriction fragment ends, followed by the calculation of a 4C-score per valid restriction fragment, taking into account if reads could be mapped to both restriction fragments ends or only one of the two [START_REF] David | HTSstation: a web application and open-access libraries for high-throughput sequencing data analysis[END_REF]. The score of the 5 restriction fragments centered on the viewpoint are set to zero. This generates the RAW 4C-seq data.

5. Smoothing and normalization of the RAW 4C-seq data. Given that RAW 4Cseq scores can be influenced, to a significant extent, by the distance between the recognition sites of the first and second restriction enzymes [START_REF] Noordermeer | The dynamic architecture of Hox gene clusters[END_REF], the base composition at its extremities, as well as differences in chromatin accessibility that could affect the digestion and ligation steps, RAW 4C-score are smoothed using a running mean across several restriction fragments. In our hands, averaging over 11 restriction fragments generates reliable 4C-seq profiles.

To compare different samples, smoothed 4C-seq profiles needs to be normalized. To do so, we consider that the sum of interaction frequencies within a region of several Megabases (e.g., a number of Topologically Associating Domains) centered on the viewpoint should be identical across different samples. We then multiply the 4C-scores by an appropriate scaling factor to obtain a sum of 1 million interaction frequency within this region.

To compare interaction frequencies between samples, a log2 ratio can be calculated between the normalized smoothed 4C-seq scores from both samples. In particular, this is done when comparing interaction frequencies between alleles, using the 4C-score from a F1 hybrid (one allele analyzed) and the one from the F1 of the reciprocal cross (second allele analyzed).

Notes

1. Different combinations of restriction enzymes can be used for high-resolution 4Cseq. In our hand, the most reliable ones for the first digestion are DpnII (and its isoschizomer MboI) and NlaIII. For optimal ligation and circularization of the targeted viewpoint, the distance between restriction sites should ideally be as follows: (a) at least 700 bp between the two consecutive recognition sites of the first restriction enzyme; (b) at least 250 bp between the recognition sites of the first and the second restriction enzyme. Some deviation from this guidance can be tolerated, especially if one of the restriction sites of the first cutter should coincide with a SNP for allelic 4Cseq. However, PCR amplifications of 4C material might be less reliable if distances become too small. [START_REF] Keane | Mouse genomic variation and its effect on phenotypes and gene regulation[END_REF] catalogs genetic variations among numerous strains of mouse. Although the official release dates from 2015 (REL-1505), two more recent releases can be downloaded from ftp://ftp-mouse.sanger.ac.uk/. The list of available strains can be queried as follows: bcftools query -l ftp://ftpmouse.sanger.ac.uk/[DIRECTORY]/[FILE.vcf.gz]. The variants of some strains are only available in the more recent releases. This is the case for the JF1_MsJ strain.

The mouse Genomes Project

3. The quality of 4C-seq primers needs to be evaluated. To validate the amplification robustness, we start by testing 4C primers lacking the Illumina adapters by running two replicate PCR reactions on agarose gel. Reliable amplifications are characterized by (1) a reproducible and uniform long smear ranging from 100 bp to 1-2 kb and (2) potentially a few discrete bands, including one that may correspond to the undigested fragment whose length can be predicted. Primers passing this first quality check are ordered with Illumina adapters and tested likewise. Primers will be rejected if the smear is too faint (hard to see on agarose gel), not sufficiently wide-ranging (for instance if centered on 200-300 bp), if many discrete bands are observed (more than 3) or if strong primer dimers are observed (not to be confused with unincorporated primers, which for sequencing material will be removed using SPRI-beads).

4.

Illumina adapters containing primers should be PAGE-purified to eliminate shorter DNA molecules aborted during oligonucleotide synthesis. To avoid the use of very long primers, a two-steps PCR strategy can be considered, whereby the first step is done with viewpoint-specific primers with a non-complementary 5'-end tail while the second amplification step relies on primers annealing to this tail and containing the Illumina adapters in 3' [START_REF] Krijger | 4C-seq from beginning to end: A detailed protocol for sample preparation and data analysis[END_REF].

5.

This protocol has been optimized for mammalian cells grown in tissue culture. When working with organism with lower DNA content per cells, the amount of starting material should be adjusted. Importantly, cell lysis may need to be optimized for each cell types, especially if starting from tissue samples. Different lysis conditions are proposed in [START_REF] Matelot | Determination of High-Resolution 3D Chromatin Organization Using Circular Chromosome Conformation Capture (4C-seq)[END_REF] [START_REF] Matelot | Determination of High-Resolution 3D Chromatin Organization Using Circular Chromosome Conformation Capture (4C-seq)[END_REF].

6.

Aggregates of nuclei might be visible at this stage. To limit their formation, reactions need to be agitated, especially during the overnight digestion with the first restriction enzyme. Using a Thermomixer and shaking at 900 rpm limits formation of aggregates.

Pipetting up-and-down the reaction with a p200 may help to dissociate the aggregates.

The presence of large aggregates is often associated with poor cell lysis and tends to inhibit digestion efficiency.

7.

DTT is unstable at ambient temperature and oxidized by oxygen. Stock solution should be aliquoted and kept frozen. When oxidized, DTT exhibits a strong absorbance at 280 nm, which alters precise DNA quantifications based on UVabsorbance. For this reason, we strongly advise to rely on fluorescence-based dyeincorporation assays for DNA quantification, especially since the DTT contained in ligation buffer will not be completely removed by the clean-up steps after the second ligation reaction. captured 3D-partners, as only one partner can be captured per viewpoint containing restriction fragment. Because the chromatin digestion is never 100% efficient, and because local 3D-partners will be captured several times, the effective number of unique 3D-partners identified should generally be expected to be smaller. The amount of 4C material used for PCR amplification can be adjusted when working with smaller or bigger genomes.

9.

Alternatively, 4C-seq PCR material can be purified with Qiagen PCR purification columns instead of SPRI magnetic beads (e.g., AMPure XP beads). Purification using a 0.8X SPRI beads ratio has the advantage that larger fragments can be removed (~200 bp for SPRI beads versus ~120 bp for PCR purification columns), which more efficiently removes unincorporated 4C-seq primers and primer dimers which tend to be incompletely removed by silica columns.

10. The multiplexing of 4C-seq PCR material originating from multiple viewpoints or samples can be done at two levels before sequencing:

-Multiple 4C-seq PCRs from different viewpoints can be mixes in equimolar quantities and sequenced together, as the first sequenced base pairs will correspond to the primer used for amplification. The demultiplexing is done by comparing the 5'-end of the read with the sequences of the different 4C-seq primers used.

-Multiple 4C-seq PCRs from different samples but with the same primer-set can also be multiplexed as the adapters used here are compatible with Illumina singleindexed sequencing. To do so, an Illumina TruSeq index is included in the P7containing primer and read using Illumina i7 Index Primer. In this case, demultiplexing is done using the dedicated Illumina software.

Up to around 20 4C-seq PCR samples, either from different viewpoints or carrying different indexes, can be pooled and sequenced together on a NextSeq 500/550 high output flow cell (aiming for 15-20 million reads per 4C-seq PCR). This number can be adjusted if using a different flow cell / sequencing platform.

11.

To ensure sufficient complexity in the first bases during the sequencing, 20% Phi-X balancer DNA should be added prior to loading on the Illumina flow cell. prevents the inclusion of the binding site for the reverse primer in the circularized ligation product (pink circle) [START_REF] Splinter | The inactive X chromosome adopts a unique three-dimensional conformation that is dependent on Xist RNA[END_REF][START_REF] Holwerda | Allelic exclusion of the immunoglobulin heavy chain locus is independent of its nuclear localization in mature B cells[END_REF]. In (B) the SNP perturbs the recognition site for the first restriction enzyme (RE1), which results in a circularized ligation product that contains the 'undigested' restriction fragment directly downstream of the viewpoint (pink circle) [START_REF] Llères | CTCF modulates allele-specific sub-TAD organization and imprinted gene activity at the mouse Dlk1-Dio3 and Igf2-H19 domains[END_REF]. In (C) a SNP is located in-between the forward 4C primer and the first cutter [START_REF] Mugford | Evidence for local regulatory control of escape from imprinted X chromosome inactivation[END_REF], which will be included in the sequencing result and thus allows the distinction between the alleles, and in (D) the SNP is located in-between the reverse 4C primer and the second cutter, which will be included in the sequencing result if paired-end sequencing is used [START_REF] Holwerda | Allelic exclusion of the immunoglobulin heavy chain locus is independent of its nuclear localization in mature B cells[END_REF][START_REF] Williams | fourSig: a method for determining chromosomal interactions in 4C-Seq data[END_REF][START_REF] De | The pluripotent genome in three dimensions is shaped around pluripotency factors[END_REF][START_REF] Dixon | Chromatin architecture reorganization during stem cell differentiation[END_REF]. Strategies (A) and (B) generate 100% allele-specific data (blue circles only), but only one of the two alleles is analyzed at a time. The other allele is analyzed using another set of primers or in the reciprocal cross if working with laboratory strains. Strategies (C) and (D) are more flexible, but the distinction between the alleles relies on the sequencing of a single polymorphism, which is sensitive to sequencing errors. For strategy (C), the SNP needs to be located in close proximity to RE1, as sequencing starts from the 5'-side of the 4C primer facing RE1. This constraint is alleviated in strategy (D), but this requires paired-end

Figures legends

4 .

 4 Immediately quench the crosslink reaction by adding 1.8 ml of a cold 1 M Glycine solution. Mix by inversion and keep the tube on ice. 5. Centrifuge (5 minutes, 300 g, 4 °C) and discard the supernatant. 6. Resuspend the cell pellet in 10 ml cold PBS. Centrifuge (5 minutes, 300 g, 4 °C) and remove most of the supernatant. 7. Resuspend the cell pellet in the remaining liquid and transfer to a 1.5 ml Safe-Lock tube.

3 .

 3 Reverse the crosslink by incubating for 1 hour and 30 minutes at 65°C while shaking in a Thermomixer (750 rpm).

4 .

 4 Digest proteins by adding 5 μl of 1 M Tris-HCl pH7.5, 2 μl of 0.5 M EDTA and 2 μl of 20 mg/ml Proteinase K. Incubate for 1 hour at 45°C while shaking in a Thermomixer (750 rpm).

6 .

 6 Transfer aqueous phase to a new 1.5 ml Safe-Lock tube. 7. Precipitate the DNA by adding 0.1 volume of 2 M NaAc pH 5.6 and 2.5 volume of 100% Ethanol. 1 µl of glycogen can be used to increase DNA precipitation efficiency. Centrifuge (30 minutes, 16,000 g, 4°C) and carefully discard the supernatant.

3 .

 3 Digest proteins of both tubes by adding 50 μl of 20 mg/ml Proteinase K and 30 μl of 20% SDS. Incubate 4 h at 45°C in a Thermomixer (900 rpm).4. Add 65 μl 5 M NaCl and incubate 1 hour at 65°C in a Thermomixer (900 rpm).5. Add 15 μl of 10 mg/ml RNAse A and incubate for 1 hour at 37°C in aThermomixer (900 rpm).6. Combine both tubes in a 15 ml Flacon tube and add 4.5 ml water to dilute out the salt of the ligation buffer that would otherwise co-precipitate in later DNA purification steps.7. Add 6 ml Phenol-chloroform-IAA to the tube, shake vigorously and centrifuge for (15 minutes, 2,200 g, room temperature). Phenol is toxic: work in a chemical fume hood with appropriate protection. 8. Transfer the aqueous phase to a 50 ml conical tube.

9 .

 9 Precipitate DNA by adding 0.1 volume of 2 M NaAc pH 5.6 and 2.5 volumes of 100% Ethanol. 5 µL Glycogen can be added to increase DNA precipitation. Mix gently and store at -80 °C for at least 1 hour. 10. Centrifuge (45 minutes, >2,200 g, 4°C) and carefully discard the supernatant. 11. Wash the precipitate with 10 ml cold 70 % Ethanol and centrifuge (10 minutes, >2,200 g, 4°C). Discard the supernatant carefully.

3 .

 3 Transfer DNA to a 50 ml conical tube. 10. Add 1 ml of Qiagen PB loading buffer to both conical tubes containing 200 µL resuspended DNA. Mix by pipetting up and down a few times. 11. Load resuspended DNA -PB buffer mix on a total of 4 QIAquick Spin Columns per sample (2 columns per conical tube).

1 .

 1 Assemble in four individual PCR tubes, with each containing a different amount of 4C material: a. 5 µL 10 X Expand Long Buffer 1 b. 1 µL 10 mM dNTPs c. 0.5 µL 100 µM forward 4C primer d. 0.5 µL 100 µM reverse 4C primer e. 4C material (25 ng, 50 ng, 75 ng or 100 ng) f. 0.75 µL (3.75U) Expand Long Enzyme mix g. Water to 50 µL 2. Run in a thermocycler as follows a. 94°C, 2 minutes b. [94°C, 15 seconds; 55-60°C, 1 minute; 68°C, 3 minutes] x 30 cycles c. 68°C, 7 minutes 3. Load 15 µL of the PCR reaction on a 1.5% Agarose gel (see Note 3)

10 .

 10 Resuspend the beads in 40 µL 10 mM Tris pH 7.5 by pipetting up-and-down at least 10 times. Incubate 5 minutes at room temperature. 11. Briefly spin the tube and place the tube on the magnetic rack. Once the solution is clear, transfer most of the solution in a new 1.5 ml microcentrifuge tube.

8 .

 8 It is critical to perform multiple 4C-PCR reactions in parallel, to ensure enough complexity before sequencing. For mammals, we amplify at least 1 µg 4C material, usually by doing 16 PCR reactions in parallel, each containing 75 ng of 4C material (16 x 75 = 1.2 µg). 1.2 µg of DNA corresponds to 4 x 10 5 copies of haploid mouse/human genome thereby allowing a theoretical maximum detection of 4 x 10 5

Figure 1 -

 1 Figure 1 -Workflow for in-nucleus Circular Chromosome Conformation Capture

Figure 2 -

 2 Figure 2 -Different strategies for allele-specific 4C-seq studies. To detect allele-

  sequencing. RE1/2: restriction site of the first and second restriction enzyme. The arrows represent the 4C primers amplifying a single (A) or both alleles (B,C,D). The yellow hexagons indicate the location of the SNP.

Figure 3 -

 3 Figure 3 -Design of the 4C-seq primers with Illumina adapters. The primers are

Figure 4 -

 4 Figure 4 -Validation of chromatin digestion by the first restriction enzyme. The

  Figure 1

2.2 Reagents and solutions

	7. DNA dye-incorporation assay for low concentration DNA (e.g., QuBit dsDNA
	High Sensitivity assay kit; Thermo Fisher Scientific, #Q32851)
	8. SPRI beads for size selection (e.g., AMPure XP; Beckman, # A63880)
	9. ATP (e.g., Sigma-Aldrich, #A2383)
	10. dNTPs (e.g., Sigma-Aldrich, # DNTP100-1KT)
	11. 20% Sodium Dodecyl sulfate (SDS) solution (e.g., Bio-Rad, #1610418)
	12. 20 mg/ml Proteinase K (e.g., NEB, # P8107S)
	13. 10 mg/ml RNase A (e.g., Thermo Fisher Scientific, #EN0531)
	14. Phenol-chloroform-IAA solution (e.g., Thermo Fisher Scientific, #15593049)
	15. 70%, 80% and 100% Ethanol solutions
	16. 10X Ligation Buffer (e.g., Promega, #C1263)
	17. 20 mg/ml BSA (e.g., NEB, #B9000S)
	18. 1 M glycine solution
	19. 5 M NaCl solution
	20. 1 M Tris-HCl, pH7.5 and 1 M Tris-HCl, pH 8 solutions
	21. 0.5 M EDTA, pH 8 solution
	22. 1 M MgCl2 solution
	1. PBS, pH 7.2 23. 2 M sodium acetate (NaAc), pH 5.6 solution
	2. Fetal Bovine Serum (FBS) 24. 1 M DTT (aliquot and store et -20°C)
	3. 37% w/v formaldehyde (e.g., Sigma-Aldrich, #47608) 25. 10X Ligation buffer (for second ligation): 660 mM Tris pH 7.5, 50 mM MgCl2,
	4. NP-40 (e.g., Sigma-Aldrich, #74385) 50 mM DTT, 10 mM ATP. Aliquot and store at -20°C for up to 6 months.
	5. Protease inhibitors (e.g., Complete protease inhibitor cocktail, EDTA free
	EASYpack; Roche, #11873580001)
	6. DNA dye-incorporation assay for high concentration DNA (e.g., QuBit dsDNA
	broad range assay kit; Thermo Fisher Scientific, #Q32850)

10. Pipettes, filter tips

2.3 Buffers, prepared within minutes before use

	2. Crosslink buffer: 2% formaldehyde solution in Cell isolation buffer.
	Formaldehyde is toxic: handling should be done in a chemical fume hood with
	appropriate protection.
	3. Cold cell lysis buffer: 10 mM NaCl, 10 mM Tris-HCl pH 8, 0.2% v/v NP-40 and
	1X Complete protease inhibitor cocktail EDTA free. Keep on ice.
	4. 0.5% SDS Solution: dilute 1 volume of 20% SDS in 39 volumes of water.
	2.4 Enzymes
	1. High concentrated restriction enzymes: DpnII (e.g., NEB #R0543M), MboI (e.g.,
	NEB #R0147M) or NlaIII (e.g., NEB #R0125L) (see Note 1)
	2. High concentrated T4 DNA ligase (e.g., Promega, M1794)
	3. Expand Long Template PCR System (Roche, # 11759060001).
	1. Cell isolation buffer: 10 % filtered FBS in PBS

Acknowledgments

We acknowledge funding from the Agence Nationale de la Recherche (project "IMP-REGULOME" -ANR-18-CE12-0022-02) to D. N.