Laure Mercier De Lépinay

Benjamin Pigeau

Benjamin Besga

Pascal Vincent

Philippe Poncharal

Olivier Arcizet

Universal Vectorial and Ultrasensitive Nanomechanical Force Field Sensor

Miniaturization of force probes into nanomechanical oscillators enables ultrasensitive investigations of forces on dimensions smaller than their characteristic length scale. Meanwhile it also unravels the force field vectorial character and how its topology impacts the measurement.

Here we expose an ultrasensitive method to image 2D vectorial force fields by optomechanically following the bidimensional Brownian motion of a singly clamped nanowire. This novel approach relies on angular and spectral tomography of its quasi frequency-degenerated transverse mechanical polarizations: immersing the nanoresonator in a vectorial force field does not only shift its eigenfrequencies but also rotate eigenmodes orientation as a nano-compass. This universal method is employed to map a tunable electrostatic force field whose spatial gradients can even take precedence over the intrinsic nanowire properties. Enabling vectorial force fields imaging with demonstrated sensitivities of attonewton variations over the nanoprobe Brownian trajectory will have strong impact on scientific exploration at the nanoscale.

Introduction-Nanosciences were revolutionized by the invention of atomic force microscopy [START_REF] Binnig | Atomic force microscope[END_REF] which enabled first perceptions of the nanoworld, by measuring proximity forces exerted by a sample on a micromechanical oscillator. The subsequent emergence of nanomechanical oscillators and evolutions in readout techniques [START_REF] Cleland | Foundations of Nanomechanics[END_REF][START_REF] Ekinci | Nanoelectromechanical systems[END_REF][START_REF] Schwab | Putting mechanics into quantum mechanics[END_REF] lead to impressive improvements in force sensitivity [START_REF] Moser | Ultrasensitive force detection with a nanotube mechanical resonator[END_REF], enabling detection of collective spin dynamics [START_REF] Degen | Nanoscale magnetic resonance imaging[END_REF][START_REF] Nichol | Nanomechanical detection of nuclear magnetic resonance using a silicon nanowire oscillator[END_REF][START_REF] Peddibhotla | Harnessing nuclear spin polarization fluctuations in a semiconductor nanowire[END_REF], single electron spin [START_REF] Rugar | Single spin detection by magnetic resonance force microscopy[END_REF], mass sensing of atoms [START_REF] Jensen | An atomic-resolution nanomechanical mass sensor[END_REF][START_REF] Sage | Neutral particle mass spectrometry with nanomechanical systems[END_REF] or inertial sensing [START_REF] Hanay | Inertial imaging with nanomechanical systems[END_REF]. Attractive perspectives arise too when nanoresonators are hybridized to single quantum systems, such as molecular magnets [START_REF] Ganzhorn | Strong spin-phonon coupling between a single-molecule magnet and a carbon nanotube nanoelectromechanical system[END_REF], spin or solid states qubits [START_REF] Arcizet | A single nitrogen-vacancy defect coupled to a nanomechanical oscillator[END_REF][START_REF] Yeo | Strain-mediated coupling in a quantum dot-mechanical oscillator hybrid system[END_REF][START_REF] Montinaro | Quantum dot opto-mechanics in a fully self-assembled nanowire[END_REF][START_REF] Pigeau | Observation of a phononic mollow triplet in a multimode hybrid spin-nanomechanical system[END_REF]. This reduction of the probe size naturally motivates the exploration of force fields on dimensions smaller than their characteristic length scale where a great physical richness is expected, in particular for nanoscale imaging or investigations of fundamental interactions such as proximity forces or near field couplings. In this situation, measurements are performed in presence of strong force field gradients whose vectorial structure becomes crucially relevant and should be fully accounted to describe the measurement process itself [START_REF] Gloppe | Bidimensional nano-optomechanics and topological backaction in a non-conservative radiation force field[END_REF]. However, most of existing force sensing scanning probe experiments are based on oscillators such as cantilevers, whose motion remains quasi-monodimensional, even in case of multimode operation [START_REF] Garcia | The emergence of multifrequency force microscopy[END_REF]. In this situation, only scalar quantities can be measured, giving access in general to the projection of the force field gradient along the measured eigenmodes which only provides partial information on the force field under investigation. In this article, we expose a novel ultrasensitive measurement technique going beyond scalar measurements and permitting to directly image vectorial 2D force fields. It is based on a singly clamped nanowire (NW) which can oscillate along two perpendicular transverse directions [START_REF] Gloppe | Bidimensional nano-optomechanics and topological backaction in a non-conservative radiation force field[END_REF][START_REF] Siria | Electron fluctuation induced resonance broadening in nano electromechanical systems: The origin of shear force in vacuum[END_REF][START_REF] Nichol | Displacement detection of silicon nanowires by polarization-enhanced fiber-optic interferometry[END_REF][START_REF] Gil-Santos | Nanomechanical mass sensing and stiffness spectrometry based on two-dimensional vibrations of resonant nanowires[END_REF][START_REF] Ramos | Optomechanics with silicon nanowires by harnessing confined electromagnetic modes[END_REF][START_REF] Cadeddu | Time-resolved nonlinear coupling between orthogonal flexural modes of a pristine gaas nanowire[END_REF]. Its vibrating extremity is immersed in the force field under investigation while its 2D-Brownian motion is optically detected and reconstructed in realtime. A novel measurement principle, based on the angular and spectral tomography of the NW Brownian motion gives access to the eigenfrequencies and eigenfunctions of the two quasi-degenerated fundamental vibrational modes. The latter are hybridized by the external force field gradients, which can be subsequently deduced from a direct mathematical inversion of the tomographic parameters after having understood the bidimensional coupling mechanism. The strength of the method relies on the use of quasi frequency-degenerated nanowires, whose eigenmodes are likely to be cross-coupled by the shear components of the external force field, generating large changes in eigenmode orientations as the force gradients easily take precedence over the nanowire intrinsic mechanical properties. We first introduce the principle of angular and spectral tomography (AngSTom) of the nanowire 2D thermal noise and demonstrate its capability to probe the in-plane 2D electrostatic force field gradients generated by a voltage-biased nano-tip approached in the vicinity of the nanowire extremity. The choice of using point source-like attractive central force field permits to fully explore the panorama of dual mode coupling in 2D which is at the heart of the measurement protocol. Finally, we realize a cartography of the vectorial structure of the electrostatic force field generated by the tip, which illustrates the strength of the method for future applications in vectorial force field scanning probe imaging. Similar NWs were previously employed [START_REF] Gloppe | Bidimensional nano-optomechanics and topological backaction in a non-conservative radiation force field[END_REF] to analyze the optomechanical backaction in 2D generated by a tightly focused laser field. To do so the optical force field was measured by means of a coherent pump-probe response measurement where the laser intensity was modulated to resonantly drive the NW with optical forces. The force measurements were realized at low optical powers, so that the nanowire dynamical response was only weakly perturbed. Beyond the demonstrated novelty, this methodology was lacking of universality since the vast majority of forces under investigation in nanoscience cannot be time modulated and a novel approach which does work in presence of strong force field gradients, as expected at the nanoscale, had to be developed. The strategy consists in exploiting a quasi frequencydegenerated nanowire, whose eigenmodes get dressed by the external force field under investigation. The dimensionality of the problem requires a novel analysis of mechanical mode coupling phenomenology in 2D which is at the heart of the method exposed.

Setup-Our nanoresonator is a 70 µm-long Silicon Carbide NW of 200 nm diameter suspended at the extremity of a sharp tungsten tip. Its vibrating extremity is immersed in a strongly focused laser beam used to probe its vibrations. The focusing is achieved by a 100x/0.75 microscope objective, producing a 330 nm optical waist. Its long working distance (4 mm) bestows a convenient buffer space between the NW and the first optical element and leaves a large spatial access to the vibrating extremity of the NW for approaching the sample under test. In this work, we investigate the electrostatic in-plane force field exerted on the NW by a 17 µm high platinum-coated pyramidal tip with a 100 nm radius of curvature. While the whole setup is grounded, this tip is electrically insulated from the rest of the setup and can be voltage biased up to a few volts. The tip position can be piezo-scanned in the vicinity of the NW extremity (see Fig. 1b) to realize force field cartography. The experiment is carried out in a vacuum chamber operated at a static pressure below 10 -2 mbar to suppress acoustic damping. Its temperature is stabilized within 0.01 • C, which limits thermal drifts below 2 nm/h. Bidimensional optomechanical readout-An optomechanical readout is employed to fully characterize the vibration properties of the nanowire in 2D and to access the spatial structure of its Brownian motion. To do so, the light backscattered on the NW extremity is collected and focused onto an amplified split photodiode which provides the sum/difference (V ⊕ /V ) of each quadrant output voltages which permits a precise micro-positioning of the NW with respect to the laser beam in the waist area, see Fig. 1c. The vibrations δr(t) of the NW extremity in the horizontal oscillating plane around its rest position r 0 are converted into temporal fluctuations in the photodetector outputs as V ,⊕ (r 0 + δr(t)). Thus both signals convey a projective measurement of the NW deflection δr ,⊕ (t) ≡ δr(t)•e ,⊕ along a measurement vector defined by e ,⊕ ≡ ∇V ,⊕ |∇V ,⊕ | evaluated at the nanoresonator rest position r 0 . By selecting a working position where both spatial gradients are at least non-collinear and ideally perpendicular (highlighted in Fig. 1c by ), a fully 2D vectorial readout of the nanowire deflection δr(t) in realtime can be achieved. Moreover, this working point lies on the optical axis which ensures the nullity of the non-conservative components of the laser dynamical backaction [START_REF] Gloppe | Bidimensional nano-optomechanics and topological backaction in a non-conservative radiation force field[END_REF], which is further circumvented by working at sufficiently low optical powers (≈ 20 -60 µW). Once the optical working position is determined, a 2D feedback loop is activated to stabilize the optomechanical probe against thermal drifts and static deflections (see SI). Nanowire 2D thermal noise-The 2D Brownian motion trajectory reconstructed from the simultaneous acquisition of δr ,⊕ (t) is shown in Fig. 1d A slight asymmetry in the NW geometry lifts the frequency degeneracy between the two mechanical polarizations to Ω 2 -Ω 1 ≈ 2π × 400 Hz whose orientations e 1,2 can be visualized in the 2D representation under resonant actuation of the NW, see Fig. 1d. To fully exploit the NW sensing capacity, it is necessary to analyze its 2D thermal noise and detect modifications of effective temperatures, eigenfrequencies and eigenmode orientations. To do so the recorded 2D NW thermal noise trajectory δr(t) is numerically projected along an arbitrary unitary measurement vector e β , providing the projected displacement δr β (t) ≡ δr(t) • e β , whose spectral density is computed and evaluated at the eigenfrequencies yielding

S δr β [Ω 1,2 ].
Repeating this sequence for varying measurement angles allows to reconstruct the angular and spectral tomography of the nanowire thermal noise in 2D, as shown in Fig. 1g. This representation permits an identification of the eigenmode orientations e 1,2 as measurement directions which maximize the noise peak heights and an estimation of the NW effective masses and temperatures, see SI.

2D force field imaging-The modification of the nanowire 2D thermal noise in presence of an external force field gradient is at the heart of the measurement principle, which is now exposed in detail. The NW dynamics, restricted to the 2 fundamental eigenmodes, is described by a 2D point-like oscillator dynamics:

δr = -Ω 2 • δr -Γδ ṙ + F(r 0 + δr) M eff + δF th M eff , (1)
where

Ω 2 ≡ Ω 2 1 0 0 Ω 2 2
is the intrinsic restoring force matrix in the e 1 , e 2 basis, Γ the mechanical damping rate, M eff their effective mass and δF th represents the 2D Langevin force vector, independently driving each eigenmode [START_REF] Pinard | Effective mass in quantum effects of radiation pressure[END_REF] with a white force noise of spectral density S th F = 2M eff Γk B T . The use of the normal mode expansion is justified by the low ( 10%) damping asymmetry observed on the uncoupled eigenmodes [START_REF] Schwarz | Deviation from the normal mode expansion in a coupled graphene-nanoresonator system[END_REF]. F(r 0 + δr) represents the action of the force field under investigation on the two fundamental eigenmodes, evaluated at the vibrating extremity position. In this work we investigate situations where the force field experienced by the NW (here of electrostatic origin) instantaneously follows -on mechanical timescales-the change in the NW position. Delays in the establishment of the force [START_REF] Ramos | Optomechanics with silicon nanowires by harnessing confined electromagnetic modes[END_REF][START_REF] Barois | Ohmic electromechanical dissipation in nanomechanical cantilevers[END_REF][START_REF] Nigues | Dynamical backaction cooling with free electrons[END_REF]] could be treated with the same formalism [START_REF] Gloppe | Generalized dynamical backaction in 2d nano-optomechancs[END_REF].

The force experienced by the NW can be linearized as

F(r 0 + δr) ≈ F(r 0 ) + (δr • ∇) F| r0 .
The first term represents the mean force undergone by the NW, responsible for a static deflection which is compensated by the 2D position feedback maintaining the NW extremity at the measurement location [START_REF] Rodrigues | Why do atomic force microscopy force curves still exhibit jump to contact?[END_REF]. As a side note, recording the correction of the feedback loop in both directions δr stat for varying positions of the sample under test (see Fig. 4c) enables a first static method to probe the 2D in-plane vectorial force field according to F = M eff Ω 2 • δr stat . A static force larger than M eff Ω 2 m k B T ≈ 1 pN will induce a nanowire deflection exceeding the spatial spreading of its thermal noise. As described below, a dynamical measurement of the NW Brownian motion around its rest position allows for a much more sensitive and robust sensing.

The second linear term in the force field expansion is responsible for a modification of the restoring forces experienced by the nanowire, which affects its dynamics and the spatial architecture of its Brownian motion. In Fourier space, the NW deflection can be expressed as δr

[Ω] = χ[Ω] • δF th [Ω]
where we have introduced the mechanical susceptibility matrix, whose inverse reads in the unperturbed basis e 1,2 :

χ -1 [Ω] = χ -1 1 [Ω] -g 11 -g 21 -g 12 χ -1 2 [Ω] -g 22 (2)
where χ -1 1,2 ≡ M (Ω 2 1,2 -Ω 2 -iΩΓ) are the unperturbed mechanical susceptibilities. Due to the immersion in the ex- ternal force field, the NW dynamics now depends on the 4 gradients of the 2D force field:

g ij (r 0 ) ≡ 1 M eff ∂ i F j | r0
whose shear components (i = j) control the crosscoupling between eigenmodes. The susceptibility matrix is diagonalized into new eigenmodes of vibration, labelled with ± indices with eigenfrequencies:

Ω 2 ± ≡ Ω 2 1 + Ω 2 2 2 ± 1 2 Ω 2 2 -Ω 2 1 2 + 4g 12 g 21 (3)
and unitary eigenvectors:

e -≡ 1 √ g 2 12 +∆Ω 2 ⊥ ∆Ω 2 ⊥ g12 and e + ≡ 1 √ g 2 21 +∆Ω 2 ⊥ -g12 ∆Ω 2 ⊥
where we have introduced

Ω 2 i ≡ Ω 2 i -g ii and ∆Ω 2 ⊥ ≡ Ω 2 2 -Ω 2 -= Ω 2 + -Ω 2 1
. The impact of the force field gradients on the NW thermal noise is illustrated in Fig. 2 in a simplified picture, while the full calculation of the 2D thermal noise spectra is carried out in the SI. The force field divergence g 11 + g 22 is responsible for a global frequency shift, while Poisson contributions g 11 -g 22 modify the frequency splitting. On the contrary, shear components g i =j are responsible for rotation of eigenmode orientations. In a conservative force field, verifying g 12 = g 21 , both eigenmodes are equally rotated, which preserves their initial orthogonality since e -• e + ∝ rot(F) • e y . The situation is qualitatively different in a non-conservative force field which breaks eigenmode orthogonality and increases their effective temperatures, see SI for full calculation. This last observation implies that a 2D vectorial readout is essential to fully determine a 2D force field. Indeed a measurement based on a single projective readout is unable to distinguish between eigenmode rotation and effective temperature increase. As can be mathematically derived, see SI, it is sufficient to determine the eigenfrequencies (Ω ± ) and eigenvectors(e ± ) of both fundamental eigenmodes to fully determine the 4 components of the force field gradients. Those are precisely the quantities that can be extracted from the spectral and angular tomography method introduced above, see SI.

Probing electrostatic force field gradients-To illustrate the strength of the method and demonstrate its validity with a well controlled force field, we employ the nanowire to probe the electrostatic force field generated by the metallic tip. The latter is piezo-positioned at a location where large shear components are generated to enhance eigenmodes hybridization. Thermal noise spectra obtained for varying bias voltage V applied on the tip are shown in Fig. 3a while the eigenfrequencies and eigenvectors orientations (e ± • e x ) are reported in Fig. 3b and3c respectively. Both eigenmodes rotate while preserving their initial orthogonality which is a direct experimental verification of the conservative character of the electrostatic force field. Strikingly, a rotation of almost 60 • can be achieved at large biases, indicating that the external force field prevails over the intrinsic NW properties to determine its eigenmode structure. The force field gradients are subsequently extracted from eigenfrequencies and eigenmode orientations and this for each bias voltage, permitting to verify their quadratic scaling according to

g ij = g 0 ij + α ij (V -V 0 ) 2 ,
as expected for electrostatic actuation, see Fig. 3. The very good agreement observed all across the bias range demonstrates the validity of the method even in extremely large force field gradients. The corresponding interpolated force field (M eff g ij δr i ) evaluated for a bias voltage of V -V 0 = 300 mV is shown in Fig. 3d. The extreme sensitivity of the NW is revealed: the relative variations of the force field magnitude over the spatial spreading of its thermal noise are only of 30 aN. The intrinsic low/high frequency eigenmodes Scanning probe measurements of electrostatic force field gradients. a Map of the force field gradients obtained for a bias voltage of -2 V while scanning the tip ≈ 500 nm below the nanowire extremity in a 6 × 6 µm 2 area (150 nm grid). b Map of relative eigenmode rotation (δθ-) illustrating the nano-compass mechanism: the e1/eigenmode is pointing towards the softer electrostatic gradient, which is generally oriented towards the nanotip. (1 µm scalebar). c Map of the 2D force field computed by integrating the measured 2D gradients for 1 V bias, using an integration constant (force vector of (-37, +4) fN) which cancels the force field just above the tip, see SI. The maximum static deflection measured was of ≈ 4 nm, see SI, which is small compared to the characteristic length over which the force field varies. d Comparison of static (blue) and dynamical (red) measurements of the force field divergence (∂xFx + ∂zFz). e Vertical cartography of the force field gradients (divergence and shear components) measured for a 1 V bias. Inset: the minimum in-plane divergence of the force field is well adjusted (red line) with a y -3 power law, see SI.

(1/2) are dressed by the external force field, defining the new set of eigenmodes (-/+) which tend to align their orientation with the eigenaxes of the external force field gradient matrix. The eigenvectors e ± ∞ of the g ij matrix are also shown and they indeed match the eigenmodes' orientations measured at large electrostatic biases. This representation also permits understanding the change of curvature visible in the frequency shifts (see SI). For example at low bias voltages ( 300 mV) the low frequency mode (1/-) still oscillates along its intrinsic orientation e 1 and experiences a trapping force field pointing towards the center, which results in a frequency increase. However at larger bias, it rotates and finally gets aligned with e - ∞ , where it experiences an anti-trapping force field, which lowers its oscillation frequency. This dynamical electrostatic force measurement [START_REF] Martin | High resolution capacitance measurement and potentiometry by force microscopy[END_REF][START_REF] Terris | Contact electrification using force microscopy[END_REF] also permits to access the local electrostatic zero, V 0 ≈ 1 V here, whose spatial variations play an essential role in surface quality analysis and proximity force measurements [START_REF] Siria | Electron fluctuation induced resonance broadening in nano electromechanical systems: The origin of shear force in vacuum[END_REF][START_REF] Behunin | Modeling electrostatic patch effects in casimir force measurements[END_REF]. Qualitatively, the longitudinal force field gradient required to shift eigenfrequencies by their mechanical linewidths amounts to ∂ 1 F 1 = k 1 /Q ≈ 10 -7 N/m, while the shear force field gradients required to rotate eigenmodes by 45

• amounts to ∂ 2 F 1 = k 1 (Ω 2 -Ω 1 )/Ω 1 ≈ 2 × 10 -7
N/m for a 0.5% intrinsic frequency splitting and a quality factor of 1000. These expressions underline the importance to work with soft NWs, presenting large aspect ratios, large quality factors and weak intrinsic asymmetries to achieve large sensitivities to both longitudinal and transverse force field gradients.

Cartography of electrostatic force field gradients-The force field gradients generated by the nanotip can be spatially mapped by piezo-scanning its position with respect to the nanowire extremity, see Fig. 4. The electrostatic tip was horizontally scanned over 40 × 40 positions in a 6 × 6 µm 2 area at a vertical distance of ≈ 500 nm below the NW extremity. In each position, the bias voltage was scanned between -2 V and 2 V while simultaneously recording the 2D thermal noise and analyzing modifications in frequency, orientations and static deflection (3 s acquisition time). The map of the low frequency mode rotation angle is shown in Fig. 4b, illustrating the nano-compass effect as this mode always rotate towards the softer direction( e - ∞ ), keeping track of the nanotip position. Reproducing in each point the analysis exposed above, the spatial map of the force field gradients ∂ i F j (r), i, j = 1, 2 can be reconstructed as shown in Fig. 4a. For our quasi point-like attractive central force field under investigation, unperturbed eigenmode orientations appear as symmetry axes in the divergence-like components (∂ i F i ) maps and as anti-symmetry axes for shear components. Shear component maps are identical (within experimental uncertainties) as expected from the conservative nature of the electrostatic actuation (∂ 1 F 2 = ∂ 2 F 1 ). A tilt in the tip scanning plane with respect to the NW axis is responsible for the slight top/bottom asymmetry observed. A phenomenological evaluation based on a radial approximation [START_REF] Derjaguin | Direct measurements of molecular attraction of solids[END_REF] is in good agreement with our measurements (see SI). The dynamically measured 2D force field gradients can be compared to the spatial gradients measured from of the static NW deflexion map (see SI): Fig 4d compares the force field 2D divergence ∝ g 11 + g 22 obtained from both methods. While both results are in quantitative agreement, as expected dynamical measurements give by far a larger measurement quality. The NW extremity can also be scanned in the vertical plane and thanks to its large sensitivity, force field gradients can still be measured microns away from the surface, see Fig. 4e. Then the 2D force field can be reconstructed through a 2D integration of the 4 force field gradients maps, including a constant force vector chosen to cancel the force field above the tip, see SI. Its attractive pattern is shown in Fig. 4c, allowing for a detailed investigation of the electrostatic behaviour of EFM tips. Beyond demonstrating the scanning probe imaging capacity of the system, this measurement also permits exploring the great richness of strong dual mode coupling in 2D: the mode crossing/anti-crossing phenomenology [START_REF] Faust | Coherent control of a classical nanomechanical two-level system[END_REF] can be tuned by carefully positioning the tip position, see SI. This also permits engineering the 2D force probe to maximize its sensitivity to specific components of the force field gradients. As an example, its sensitivity to shear component can be maximized by electrostatically minimizing its intrinsic frequency splitting. The measurements accumulated at varying distances (both horizontally and vertically) and bias voltages have permitted verifying the linearity of our probe on a large dynamical range spanning from a few nN/m to tens of µN/m.

Conclusions-

We have demonstrated a novel approach allowing vectorial force field imaging at the nanoscale, especially suited for probing in-plane 2D force fields. By bringing an additional dimension to surface imaging techniques based on mechanical signals with demonstrated force sensitivities in the attoNewton range at room temperature -several orders of magnitude better than commercial AFM apparatus-we anticipate that our approach will quickly bring novel perspectives in scientific imaging at the nanoscale. It will contribute to shed new light on the field of probe-surface interaction [START_REF] Siria | Electron fluctuation induced resonance broadening in nano electromechanical systems: The origin of shear force in vacuum[END_REF][START_REF] Mate | Atomic-scale friction of a tungsten tip on a graphite surface[END_REF][START_REF] Karrai | Interfacial shear force microscopy[END_REF][START_REF] Kuehn | Dielectric fluctuations and the origins of noncontact friction[END_REF]. Developments of novel measurement protocols based on driven 2D trajectories will permit both increasing the measurement speed and improving the sensitivity, while simultaneous observation of several longitudinal modes families [START_REF] Garcia | The emergence of multifrequency force microscopy[END_REF] will further enrich the imaging capacity of our vectorial force probe. In particular, it could be straightforwardly employed to explore proximity forces at the nanoscale, such as Casimir forces in nano-structured samples where novel phenomenology can be expected [START_REF] Rodriguez | The casimir effect in microstructured geometries[END_REF][START_REF] Chen | Demonstration of the lateral casimir force[END_REF][START_REF] Guérout | Private communication[END_REF]. The method is also compatible with non-conservative force field imaging and will permit further exploration of fluctuation theorems in 2D. Finally, our 2D analysis can also be transposed to 3D, using 3D angular and spectral analysis with for instance trapped ions [START_REF] Bushev | Shot-noise-limited monitoring and phase locking of the motion of a single trapped ion[END_REF] or levitated nano-particles [START_REF] Li | Millikelvin cooling of an optically trapped microsphere in vacuum[END_REF][START_REF] Gieseler | Subkelvin parametric feedback cooling of a laser-trapped nanoparticle[END_REF].

24 ,FIG. 1 .

 241 FIG. 1.2D optomechanical readout. a Schematics of the experiment. The transverse vibrations of a singly-clamped SiC nanowire are optically readout using a dual photodiode amplifier recording the reflexion of a laser beam tightly focussed on the NW extremity. b CCD image of a voltage biased sharp metallic tip piezo-positioned at the NW extremity to generate the static electrostatic force field under investigation. c Reflection maps of the DC differential/sum signals V ⊕/ (r) obtained by piezo-scanning the nanowire in the laser waist. The measurement location indicated by provides orthogonal projective measurement vectors, e ,⊕ allowing a 2D reconstruction of the nanowire trajectories δr(t). d Driven trajectories obtained by resonant actuation of both fundamental eigenmodes. e Thermal noise trajectory recorded over 100 ms (100 ns sampling). f Calibrated projected displacement noise spectra S δr ,⊕ [Ω]. g Angular spectral tomography the 2D thermal noise, showing the relative variations of S δr β [Ω1,2] with the analysis angle β, giving access to the eigenmode orientations e1,2.

  . It presents a gaussian distribution spreading over ∆x = k B T /M eff Ω 2 m ≈ 5 nm, which is ≈ 2.5% of the NW diameter. Corresponding vibration noise spectra S δr ,⊕ [Ω] are shown in Fig 1f, restricted to the fundamental longitudinal mode family oscillating around 77 kHz with quality factors Q ≈ 1000 at 5 × 10 -5 mbar and effective masses M eff = 600 fg (stiffness k i ≈ 150 µN/m).

FIG. 2 .

 2 FIG. 2. Measurement principle of 2D force field gradients. Any 2D force field gradients can be fully expanded on the geometric basis represented in a. Their respective impact on thermal noise spectra (measured along one single measurement vector) and on ASTom representations are shown in b and c. The knowledge of the modified eigenvectors' orientations and eigenfrequencies allows a complete derivation of the local external force field gradients.

FIG. 3 .

 3 FIG. 3.Mapping electrostatic force field gradients A voltage bias is applied on the electrostatic tip to tune the external force field gradient experienced by a ≈ 40 kHz NW for a nanotip position chosen to generate large cross-coupling between eigenmodes. a Thermal noise spectra S δr ,⊕ [Ω] obtained for varying voltage bias. b,c Deduced frequency shifts and eigenmode rotations whose knowledge permits a full derivation of the force field gradients ∂iFj. Solid lines are fits according to the theory, based on a quadratic dependence of the electrostatic force field, using M eff gij = {-146, 100, 114, 14} nN/m/V 2 (i,j=1,2) and an offset voltage of 1 V, see text and SI. d Linearly extrapolated cartography of the local force field (F(r) -F(r0)) drawn over the spatial extension of the NW thermal noise for a 300 mV bias. At large bias, the NW eigenvectors get aligned with the eigenvectors e ± ∞ of the electrostatic force field gradient matrix gij.

  FIG. 4.Scanning probe measurements of electrostatic force field gradients. a Map of the force field gradients obtained for a bias voltage of -2 V while scanning the tip ≈ 500 nm below the nanowire extremity in a 6 × 6 µm 2 area (150 nm grid). b Map of relative eigenmode rotation (δθ-) illustrating the nano-compass mechanism: the e1/eigenmode is pointing towards the softer electrostatic gradient, which is generally oriented towards the nanotip. (1 µm scalebar). c Map of the 2D force field computed by integrating the measured 2D gradients for 1 V bias, using an integration constant (force vector of (-37, +4) fN) which cancels the force field just above the tip, see SI. The maximum static deflection measured was of ≈ 4 nm, see SI, which is small compared to the characteristic length over which the force field varies. d Comparison of static (blue) and dynamical (red) measurements of the force field divergence (∂xFx + ∂zFz). e Vertical cartography of the force field gradients (divergence and shear components) measured for a 1 V bias. Inset: the minimum in-plane divergence of the force field is well adjusted (red line) with a y -3 power law, see SI.

Acknowledgements-We warmly thank the PNEC group at ILM, A. Gloppe, B. Canals, J. Chevrier, S. Reynaud, A. Lambrecht, J.P. Poizat, G. Bachelier, J. Jarreau, C. Hoarau, E. Eyraud and D. Lepoittevin, for theoretical, experimental and technical assistance. This project is supported by the ANR (FOCUS-13-BS10-0012), the ERC Starting Grant StG-2012-HQ-NOM and Lanef (CryOptics).


Jonathan Viquerat

Nikolai Schmitt

Claire Scheid

Simulating 3D periodic structures at oblique incidences with discontinuous Galerkin time-domain methods: theoretical and practical considerations

Keywords: computational electromagnetics, time-domain Maxwell equations, discontinuous Galerkin method, periodic structures, oblique incidence sources, nanophotonics

come

Introduction

A considerable amount of structures studied numerically in the nanophotonics field can be modeled using Periodic Boundary Conditions (PBC). Thus, only a unit cell of the computational domain needs to be considered, the PBC mimicking an infinite repetition of the pattern in space. Periodic structures in nanooptics cover a vast range of applications, such as color filters [START_REF] Park | Structural color filters enabled by a dielectric metasurface incorporating hydrogenated amorphous silicon nanodisks[END_REF], reflectarrays for wireless communications [START_REF] Zou | Dielectric resonator nanoantennas at visible frequencies[END_REF], or sensing devices [START_REF] Lee | Surface plasmon resonance sensing: Periodic metallic nanostructures for high-sensitivity biosensing applications[END_REF]. In many situations, a wide-band plane wave analysis at normal incidence provides enough information about the structure. However, some applications also require the knowledge of the device response for a range of incidence angles.

For a better understanding, let us suppose a 2D rectangular domain periodic in direction x and of lateral size δ x . Then, suppose a plane wave propagating in the ẑ direction (i.e. at normal incidence) in the domain. In this case, the PBC for any field U can be written as:

U(x, z, t) = U(x + δ x , z, t).
This situation is represented in figure 1(a), where the two matching periodic points (green dots) can update their field informations using each other's data at the same time-step. In the case of oblique incidence (figure 1(b)), the update of the left green dot would require information of the right green dot, which is located in an area that has not yet been hit by the incident field. In other words, this means that the field update at the left point requires the knowledge of the fields at the right point at a future timestep: and is infinite in the vertical direction. Here the notation PMLs at top and bottom already refers to the artificial numerical truncation of the domain that will be used in the sequel. The incident plane wave is oriented toward the bottom of the domain.

U(x, z, t) = U(x + δ x , z, t + δ x sin θ c 0 ), (1)
Numerically, this will be imposed at an artificial interface called the Total-Field/Scattered-Field interface (see later for details).

In normal incidence, the PBC applies naturally, and at each discretization point matching information is available from the periodic neighbor point (green dots). In oblique incidence, a causality issue appears, since the matching information for the left green dot is located in a zone that has not been illuminated by the incident field yet.

where c 0 is the speed of light in vacuum. Because of this causality issue, standard numerical methods such as Finite Difference Time Domain (FDTD) or Discontinuous Galerkin Time Domain (DGTD) applied to the governing Maxwell's equations fail to use PBC as is to simulate structures in the case of an oblique excitation. Over the years, numerous techniques have been developed in the FDTD field to deal with this modeling issue:

1. The sine-cosine method [START_REF] Kokkinos | Periodic finite-difference time-domain analysis of loaded transmission-line negative-refractive-index metamaterials[END_REF], where fields are split on two separate grids corresponding to its sine-and cosine-dependent components. In this case, a generalized PBC can be used, although this restricts the computation to single-frequency sources; 2. The sliding unit cell method [START_REF] Hotter | Broadband FDTD analysis of infinite phased arrays using periodic boundary conditions[END_REF], in which the unit cell itself is slided as the excitation progresses in the domain. This method presents several restrictions on the frequencies that can be analyzed, and several values of the field must be stored at boundaries of the domain for future updates; 3. The multiple unit cell method [START_REF] Lee | A conceptually simple method for incorporating periodic boundary conditions into the FDTD method[END_REF], where several unit cells are used at the same time in order to separate the correct field propagating in the periodic structure from the error due to the causality issue. On top of requiring field storing for future updates, this method requires two unit cells in 2D, and four in 3D, which mitigates the interest of PBC; 4. The angle-updated method [START_REF] Taflove | Computational Electrodynamics: The Finite-Difference Time-Domain Method, 3rd Edition[END_REF], in which the time-domain update scheme is rearranged in order to compute a part of the fields in the domain before starting the actual FDTD computation. This method suffers from a severe limitation in the incidence angles that can be used (θ < 45 • in 2D, and θ < 35 • in 3D); 5. The field-transform and the split-field techniques [START_REF] Roden | Time-domain analysis of periodic structures at oblique incidence: Orthogonal and nonorthogonal FDTD implementations[END_REF]: first, a field transformation is applied to the continuous Maxwell's equation that cancels out the causality problem. Then, for compatibility with the staggered grid inherent to Finite Difference, the fields are split before discretization. This method is well spread among the FDTD community, since it allows for any incidence angle below 90 • [START_REF] Taflove | Computational Electrodynamics: The Finite-Difference Time-Domain Method, 3rd Edition[END_REF]. However, its counterpart is a reduction of the maximal stable timestep in the time discretization with the incidence angle (see next section).

Although a large amount of references are available for such techniques in the FDTD framework (see [START_REF] Valuev | Implementation of the iterative finite-difference time-domain technique for simulation of periodic structures at oblique incidence[END_REF] and references therein), for the DGTD method the literature remains very shallow. To the best knowledge of the authors, the only methodological contribution on this topic is from Miller et al [START_REF] Miller | A discontinuous Galerkin time domain framework for periodic structures subject to oblique excitation[END_REF], although the method has been exploited afterwards [START_REF] Zeng | Modeling and design of near-field antennas with periodic structures[END_REF]. In [START_REF] Miller | A discontinuous Galerkin time domain framework for periodic structures subject to oblique excitation[END_REF], the authors revisit the field transform method, and then present its discretization in the DGTD framework. No analysis is provided, neither for the continuous nor for the discrete equations. Regarding domain truncation, only first-order Absorbing Boundary Conditions (ABC) are covered, and very little is said about the incident fields or the observables.

In this paper, we focus on developing the field transform technique in the context of nanophotonics periodic structures. We use a DGTD discretization framework for the discretization of the governing Maxwell's equation and provide a stability analysis. We extend the technique to incorporate the possibility of using PML's and finally, we give several numerical results from academical to more applied test cases.

In the following section, we shortly sum up the field transform technique, and the continuous system is analyzed with the method of characteristics and an energy estimate is provided. Then, its discretization with the DGTD method is covered, and a stability proof is given in this context. Along the way, insights about the transformation of sources, observables (reflectance, transmittance and diffraction efficiency), and the use of CFS-PMLs for the transformed system are given and included in the analysis. We then consider the numerical implementation and validate the latter on a textbook case for which an analytical solution is known. Finally, we consider two realistic cases: (i) a silicon-based color filter, which response is analyzed under various incidence angles, and (ii) a diffraction grating, for which we compare the diffraction efficiencies at normal and oblique incidences.

The field transform technique

Transformed equations

We denote ε 0 , µ 0 , c 0 respectively the permittivity, permeability and speed of vacuum and follow [START_REF] Miller | A discontinuous Galerkin time domain framework for periodic structures subject to oblique excitation[END_REF] for the setting of the equations. Consider a plane wave of angular frequency ω defined by a wave vector k = t (k x , k y , k z ) traveling in a background material where its wavefront speed is c. We introduce c b the relative speed of the wave in this material, such that c = c b c 0 , and denote ε b and µ b the corresponding relative electromagnetic parameters. The permittivity and permeability of the background material are thus respectively given by ε = ε b ε 0 , µ = µ b µ 0 . We define a normalized wave vector:

k = k |k| = t (κ x , κ y , κ z ) = t (sin θ cos φ, sin θ sin φ, cos θ) ,
where |k| = ω c , and (θ, φ) are defined as the classical spherical angles, see figure 2. The vector k is decomposed in k = k + k⊥ , where k⊥ is oriented in the common direction of the periodicity planes, which is here taken to be ẑ for the sake of simplicity. Then: k = t (κ x , κ y , 0) and k⊥ = t (0, 0, κ z ) .

We then consider any periodic structure with electromagnetic constants ε = ε 0 ε s and µ = µ 0 µ s (ε s and µ s being the relative parameters) embedded in the background material. Let us now recall the source-free frequency-domain Maxwell's equations in a media of electromagnetic constants (ε, µ):

jωεE(r, ω) = ∇ × H(r, ω), jωµH(r, ω) = -∇ × E(r, ω).
(

) 2
The field transform technique consists in formally considering the following transformation:

P(r, ω) = E(r, ω) exp jk • r , S(r, ω) = H(r, ω) exp jk • r . (3)
Applying transformation (3) to system (2), a few lines of calculations along with the identity ∇ × (ψV) = (∇ψ) V + ψ (∇ × V) lead to: jωεP(r, ω) = ∇ × S(r, ω) -jk × S(r, ω), jωµS(r, ω) = -∇ × P(r, ω) + jk × P(r, ω).

Formally transforming back to time-domain with the identity jω ↔ ∂ ∂t and recalling that k = k |k| yields the following system in the time dependent equations:

ε ∂P ∂t (r, t) + k c 0 c b × ∂S ∂t (r, t) = ∇ × S(r, t), µ ∂S ∂t (r, t) - k c 0 c b × ∂P ∂t (r, t) = -∇ × P(r, t).
Finally, for the sake of simplicity, we classically normalize the latter system following [START_REF] Viquerat | Simulation of electromagnetic waves propagation in nano-optics with a high-order discontinuous Galerkin time-domain method[END_REF] to obtain:

ε r ∂P ∂t (r, t) + k c b × ∂S ∂t (r, t) = ∇ × S(r, t), µ r ∂S ∂t (r, t) - k c b × ∂P ∂t (r, t) = -∇ × P(r, t). (4)
with ε r = ε s , µ r = µ s in the structure and ε r = ε b , µ r = µ b in the background material. Forming the vector U = t (P, S), one can recast system (4) under a more general form:

Q ∂U ∂t + A ∂U ∂x + B ∂U ∂y + C ∂U ∂z = 0, (5)
where A, B, C are the usual coefficient matrices (see definition in [START_REF] Shang | Characteristic based methods for the time-domain Maxwell equations[END_REF]), and Q is defined as:

Table 1: Characteristic speeds for system [START_REF] Hotter | Broadband FDTD analysis of infinite phased arrays using periodic boundary conditions[END_REF], with c = c b cr .

λ i Multiplicity κx + c 2 -κ 2 y c 2 -κ 2 x -κ 2 y c b 2 κx -c 2 -κ 2 y c 2 -κ 2 x -κ 2 y c b 2 κy + c 2 -κ 2 x c 2 -κ 2 x -κ 2 y c b 2 κy -c 2 -κ 2 x c 2 -κ 2 x -κ 2 y c b 2 c b c 2 -κ 2 x -κ 2 y 2 -c b c 2 -κ 2 x -κ 2 y 2 0 6 Q = 1 c b         c b ε r 0 0 0 0 κ y 0 c b ε r 0 0 0 -κ x 0 0 c b ε r -κ y κ x 0 0 0 -κ y c b µ r 0 0 0 0 κ x 0 c b µ r 0 κ y -κ x 0 0 0 c b µ r         .
Note that system (4) both depends (i) on the local material, through ε r and µ r , and (ii) on the background material, through its relative speed of light c b . The background material should be that where the incident field is imposed. Its introduction in the transformation allows to tackle cases where the considered structure is not embedded in vacuum. The latter system can be analyzed on the structure using an eigenvalue study.

Characteristics analysis

For a hyperbolic system, we can follow the method of characteristics [START_REF] Shang | Characteristic based methods for the time-domain Maxwell equations[END_REF], use the characteristic variables (which are not detailed here) and conclude that they will travel along trajectories which slopes are the eigenvalues of Q -1 A, Q -1 B and Q -1 C. These eigenvalues are shown in table 1. Let us now analyze the behavior of these eigenvalues in a few basic cases. For the sake of simplicity, in this section, we restrict ourselves to non-magnetic materials (µ b = µ s = 1).

2.2.1. φ = 0, c = 1
Having φ = 0 (i.e. the wave will propagate parallel to the {y = 0} plane) yields κ x = sin θ and κ y = 0. Additionally, c = 1 indicates that we are looking at the behavior of the eigenvalues of system (4) in the background material. In this case, their expressions reduce to: For θ = 0 • , one recovers the regular eigenvalues λ = ±c b (which is true even for φ = 0). The latter expressions show that, in vacuum, the speed of some characteristic variables of the transformed system (4) are superior to that of light as soon as θ > 0 • . This fact is coherent with the reduction of the largest stable timestep with increasing θ noted in [START_REF] Miller | A discontinuous Galerkin time domain framework for periodic structures subject to oblique excitation[END_REF], and with our own theoretical predictions and numerical experiments that will be presented in next sections. As shown in figure 3, the most constraining speed in terms of CFL value is λ x + (θ), meaning that the maximal stable timestep will be expected to decrease as

λ x + ,x -(θ) = ±c b 1 ∓ sin θ , λ y + ,y -(θ) = ±c b 1 -sin 2 θ , λ z + ,z -(θ) = λ y + ,y -.
1 λ x + (θ) .
2.2.2. φ = 0, c = 1 Having c = 1 means that we are considering the solutions of system (4) in a material that is not the background material. The eigenvalues then have the following expressions:

λ x + ,x -(θ) = ±c b c ∓ sin θ , λ y + ,y -(θ) = ±c b c 2 -sin 2 θ , λ z + ,z -(θ) = λ y + ,y -,
The speed of the characteristic variables is still superior to that of the background material. However, two cases now need to be distinguished:

1. c > 1, i.e. ε r > ε b (for example, a dielectric slab embedded in vacuum). All values of θ in the range [0 • , 90 • [ are acceptable; 2. c < 1, i.e. ε r < ε b (
for example, a vacuum slab between two dielectric half-spaces). In this case,

λ x + (θ) is singular for a critical value θ c : θ c ( c ) = sin -1 c , (6)
while λ y + (θ) and λ y -(θ) become imaginary with opposite signs for θ > θ c . A plot of θ c as a function of c is given in figure 4. Following Snell's law, it should be noted that θ c is nothing else than the angle of total internal reflection, for which all the light will end up travelling along the materials interface.

A direct consequence is that system (4) is not stable for θ ≥ θ c if c < 1. This assumption will be confirmed with the energy study and verified numerically in the results section.

General case

We recover the critical angle, since we have to impose that:

c 2 -κ 2 x -κ 2 y > 0, which implies: sin 2 (θ) < c 2 .
Imposing this condition suffices to have both c 2 -κ 2 y ≥ 0 and c 2 -κ 2 x ≥ 0.

Energy principle

In this section, we investigate the energy of the system. The physical domain is denoted by Ω ⊂ R 3 , and is truncated in the ẑ direction. The boundary of Ω, Γ, thus consists of several parts: Γ ext the exterior boundary and Γ P , the periodic boundary. On Γ ext using a scattered field formulation, will be either imposed absorbing boundary conditions or a layer of PML's (see next section). The domain Ω itself consists in the background material Ω b (with physical parameters ε r = ε b and µ r = µ b ) and the structure Ω s (with physical parameters ε = ε s and µ = µ s ).

We denote by •, • D and • D respectively the L 2 scalar product and L 2 norm on a given set D.

We consider the equations on a fixed time interval [0, T ], T > 0. To summarize, we consider the following system of equations, for all t ∈ [0, T ]:

ε r ∂P ∂t (•, t) + k c b × ∂S ∂t (•, t) = ∇ × S(•, t), on Ω, µ r ∂S ∂t (•, t) - k c b × ∂P ∂t (•, t) = -∇ × P(•, t), on Ω, (7)
These equations are supplemented with initial and boundary conditions. For the latter, we impose PBC on the lateral boundaries of the domain (see figure 1). The outer boundaries are not specified right now to be able to envisage several truncation alternatives.

Remark 1. In this paper, the focus is put on stability analysis and numerical computations rather that any well posedness study, that would directly follow from classical results for Maxwell's equations.

Proposition 1. Consider that (S, P) verifies (7) and are smooth fields, i.e. (S, P) ∈ C 1 (0, T, H(curl, Ω)) 2 .

If for all t ∈ [0, T ],

E(t) := 1 2 ε r P(t, •) 2 Ω + 1 2 µ r S(t, •) 2 Ω + k c b × S(t, •), P(t, •) Ω . Then E(t) = E(0) + t 0 S × n, P Γext .
Proof. We take the L 2 -scalar product of ( 7) with test function (S, P). We obtain:

ε r ∂P ∂t , P Ω + k c b × ∂S ∂t , P Ω = ∇ × S, P Ω , µ r ∂S ∂t , S Ω - k c b × ∂P ∂t , S Ω = -∇ × P, S Ω ,
which is equivalent to:

1 2 ∂ ∂t ε r P, P Ω + k c b × ∂S ∂t , P Ω = ∇ × S, P Ω , 1 2
∂ ∂t µ r S, S Ω - k c b × ∂P ∂t , S Ω = -∇ × P, S Ω .
Performing an integration by parts, one finds:

1 2 ∂ ∂t ε r P, P Ω + k c b × ∂S ∂t , P Ω = S, ∇ × P Ω + S × n, P ∂Ω , 1 2
∂ ∂t µ r S, S Ω - k c b × ∂P ∂t , S Ω = -∇ × P, S Ω .
Let us detail the boundary term S × n, P ∂Ω :

S × n, P ∂Ω = S × n, P Γ P + S × n, P Γext .

Using the periodic boundary conditions on Γ P , we deduce that S × n, P Γ P = 0. Summing the two equations:

1 2 ∂ ∂t ε r P, P Ω + 1 2 ∂ ∂t µ r S, S Ω + k c b × ∂S ∂t , P Ω - k c b × ∂P ∂t , S Ω = S × n, P Γext
Using the properties of the mixed product, one has:

1 2 ∂ ∂t ε r P, P Ω + 1 2 ∂ ∂t µ r S, S Ω + k c b × ∂S ∂t , P Ω + k c b × S, ∂P ∂t Ω = S × n, P Γext Thus, 1 2
∂ ∂t ε r P, P Ω + 1 2 ∂ ∂t µ r S, S Ω + ∂ ∂t k c b × S, P Ω = S × n, P Γext ,
and the result follows.

We now prove that the energy E is positive under some assumption on k and apply this result in the special cases considered in this paper. Proposition 2. Under the hypotheses of proposition 1, one has for all t ∈ [0, T ]:

E(t) ≥ 1 2 ε r   1 - k c   P(t, •) 2 Ω + 1 2 µ r   1 - k c   S(t, •) 2 Ω . (8)
If one supposes that k < c , then the function t → E(t) is positive quadratic form on [0, T ]. Furthermore, if the exterior boundary conditions are chosen such that t 0 S×n, P Γext ds ≤ M , with M a positive constant, then:

E(t) ≤ E(0) + M, (9)
i.e. the positive energy is bounded.

Proof. One has:

k c b × S, P Ω ≤ k c b S Ω P Ω , ≤ k c c r S Ω P Ω ,
Using weighted Young's inequality, we obtain k

c b × S, P Ω ≤ 1 2 ε r k c P 2 Ω + 1 2 µ r k c S 2 Ω
This implies that:

k c b × S, P Ω ≥ -   1 2 ε r k c P 2 Ω + 1 2 µ r k c S 2 Ω  
Using this estimate, one has the following:

1 2 ε r   1 - k c   P 2 Ω + 1 2 µ r   1 - k c   S 2 Ω ≤ E(t)
The conclusions follow from this last inequality and proposition 1.

Remark 2. If one considers absorbing boundary conditions (first order Silver-Müller) or perfect electric condition on Γ ext , then S × n, P Γext ≤ 0. In this case, the energy is decreasing.

Remark 3. We recover the stability condition inferred in the characteristic analysis through the condition k ≤ c .

Remark 4. We would like to point out that this study naturally extends to classical dispersion models such as Drude model or even generalized dispersion laws that are commonly used in nanophotonics. Indeed, these models consist only in adding auxiliary ordinary differential equations (see [START_REF] Lanteri | Analysis of a generalized dispersive model coupled to a DGTD method with application to nanophotonics[END_REF] for details on the models) that are linearly coupled to Maxwell equations via a source current. For the sake of clarity, we do not detail this extension here, since it is straightforward. We refer the reader to the discussion on this topic in section 3.5).

We now turn to the numerical approximation. The discretization framework is based on a DGTD scheme. Since the type scheme that we use is well described in a significant number of references (among them see e.g. [START_REF] Miller | A discontinuous Galerkin time domain framework for periodic structures subject to oblique excitation[END_REF] in the context of PBC, [START_REF] Viquerat | Simulation of electromagnetic waves propagation in nano-optics with a high-order discontinuous Galerkin time-domain method[END_REF][START_REF] Lanteri | Analysis of a generalized dispersive model coupled to a DGTD method with application to nanophotonics[END_REF] in a nanophotonic context), we choose to give the reader a quick presentation of it and refer to the above mentioned references for a more thorough description.

The DGTD method

Numerical scheme

The steps of the derivation of the DGTD scheme of system (4) are well described in [START_REF] Miller | A discontinuous Galerkin time domain framework for periodic structures subject to oblique excitation[END_REF]. An important point is that the numerical flux, that is introduced through the DGTD scheme, does not need to be modified for the transformed equations. This makes the derivation straightforward.

We consider that Ω ⊂ R 3 is a bounded convex domain (in fact here, even parallelipipedic as described previously), and n the unitary outward normal to its boundary ∂Ω. Let Ω h be a discretization of Ω, relying on a quasi-uniform triangulation T h verifying T h = N i=1 T i , where N ∈ N * is the number of mesh elements, and (T i ) i∈ 1,N the set of simplices. For each cell T i , V i is the set of indices {k ∈ 1, N | T i T k is a triangular face}, and n ik is the unit normal to the face, oriented from T i to T k . Then, let us fix k > 0 and introduce V h , the subspace of L 2 (Ω) fields that are polynomial of degree at most k, on each element of the mesh. For each W h ∈ V h , we note W i , the vector of degrees of freedom on W h restricted to the element T i on the classical Lagrange nodal basis.

We use an element-wise weak formulation for (discontinuous) unknowns in V h , where communication within two cells of the mesh is recovered by introducing fluxes. The semi-discrete DG scheme reads in a matricial form: find (P h , S h ) ∈ V h × V h such that for all i ∈ {1, ..., N }, ∂ ∂t

P i S i = Q -1 M -1 i -K i × S i + l∈Vi S il (S * × n il ) K i × P i -l∈Vi S il (P * × n il ) , (10)
where M i is the local mass matrix on the element T i , K i the rigidity matrix on the element T i , S il the interface matrix. Furthermore, P * and S * represent the numerical fluxes on the interfaces [START_REF] Viquerat | Simulation of electromagnetic waves propagation in nano-optics with a high-order discontinuous Galerkin time-domain method[END_REF]:

P * = 1 Y i + Y l ({Y P} il + αn il × S il ) , S * = 1 Z i + Z l ({ZS} il -αn il × P il ) . ( 11
)
on T i ∩ T l , i ∈ {1, .., N }, l ∈ V i . The quantities Y and Z are respectively the material admittance and impedance. The parameter α is a given real in [0, 1]. When α = 0, one obtains the so-called centered flux. When α = 1, one obtains the fully upwind flux. Semi-discrete initial conditions are chosen in accordance with the continuous ones such as P h (0, .) = Π h (P(0, •)) and S h (0, •) = Π h (S(0, •)), with Π h a well-chosen projection operator such as the L 2 one.

Proposition 3. Let us suppose that (S h , P h ) ∈ C 1 ([0, T ], V h ) × C 1 ([0, T ], V h ) is a solution of system (10) for a given α ∈ [0, 1]. If one defines: E h (t) := 1 2 ε r P h (t, •) 2 Ω + 1 2 µ r S h (t, •) 2 Ω + k c b × S h (t, •), P h (t, •) , (12)
then E h is a positive quadratic form under the assumption that k < c . Furthermore,

E h (t) ≤ E h (0) + t 0 S h × n, P h Γext ds (13)
Proof. We rewrite [START_REF] Miller | A discontinuous Galerkin time domain framework for periodic structures subject to oblique excitation[END_REF] in its variational form. Then the proof is a direct combination of the arguments that appeared in Propositions 1 and 2 and classical arguments, depending on the fluxes used (see e.g. [START_REF] Lanteri | Analysis of a generalized dispersive model coupled to a DGTD method with application to nanophotonics[END_REF]). We choose not to reproduce it here and will rather focus on the added difficulty of the time discretization.

We directly conclude to the following Corollary 1. If boundary conditions are chosen such that there exists M > 0 (independent of the discretization parameters) such that ∀t ∈ [0, T ], t 0 S h × n, P h Γext ≤ M , then the semi-discrete scheme is stable.

Remark 5. As for the continuous case, the semi-discrete scheme is thus in particular stable if one chooses absorbing boundary conditions or perfect electric boundary conditions.

Regarding the time integration scheme of [START_REF] Miller | A discontinuous Galerkin time domain framework for periodic structures subject to oblique excitation[END_REF], we first propose to investigate the stability of a simple classical explicit scheme (leap-frog with centered fluxes) in order to understand the impact of the field transform technique on the CFL condition. In the numerical simulations, for accuracy reasons, we would preferably use the well-known five-stages fourth-order Low-Storage Runge-Kutta (LSRK) method described in [START_REF] Carpenter | Fourth-order 2n-storage Runge-Kutta schemes[END_REF] combined with upwind fluxes (α = 1).

We consider a given uniform subdivision of the time interval [0, T ] of J + 1 points (J ∈ N * ) and ∆t > 0 the corresponding time step. A leap-frog time integration scheme of [START_REF] Miller | A discontinuous Galerkin time domain framework for periodic structures subject to oblique excitation[END_REF] with centered fluxes (i.e. α = 0 in the expression of the fluxes ( 11)) gives the fully discrete problem: For all n ∈ {0, ..., J -1}, find (P n , S n+ 1 2 ) ∈ V h × V h such that on each element T i (i ∈ {1, ..., N }):

   ε r P n+1 i -P n i ∆t + k c b × S n+ 3 2 i -S n+ 1 2 i ∆t µ r S n+ 3 2 i -S n+ 1 2 i ∆t - k c b × P n+1 i -P n i ∆t    = -K i × S n+ 1 2 i + l∈Vi S il S n+ 1 2 * × n il K i × P n+1 i -l∈Vi S il P n+1 * × n il , (14)
We define the modified local energy at time step n + 1 2 on the element T i as follows:

E n+ 1 2 i = 1 2 ε r P n 2 Ti + µ r 1 2 S n+ 1 2 2
Ti -

1 4 ∆t ∇ × S n+ 1 2 , P n Ti + S n+ 1 2 , ∇ × P n Ti - 1 2 ∆t S n+ 1 2 + × n, P n ∂Ti + k c b × S n+ 1 2 , P n Ti , (15)
and the total modified energy:

E n+ 1 2 h = i E n+ 1 2 i .
Here, S + denotes S l on a face a il of T i (i.e. the value of S h on a il from T l ).

Proposition 4. The energy E n+ 1 2 is positive under a CFL condition given as

∆t h ≤ M c -k c b ,
with M a positive constant independent of the discretization parameters.

Proof. One makes use of Young's inequality and inverse inequalities stating that there exists C > 0 independent of the discretization parameters such that for all

U h ∈ V h , ∇ × U h Ti ≤ C h i U h Ti and U h ∂Ti ≤ C √ h i U h Ti
and Young's inequality (here h i denotes the diameter of T i ). This gives that there exists C > 0 such that:

E n+ 1 2 i ≥ 1 2 ε r P n 2 Ti + 1 2 µ r S n+ 1 2 2
Ti -

1 4 C c∆t h i ε r S n+ 1 2 2
Ti -

1 4 C c∆t h i µ r P n 2 Ti - 1 2 ε r k c P n 2 Ti - 1 2 µ r k c S n+ 1 2 2
Ti , which simplifies to:

E n+ 1 2 i ≥ ε r ( 1 2 - 1 4 C c∆t h i - 1 2 k c ) P n 2 Ti + µ r ( 1 2 - 1 4 C c∆t h i - 1 2 k c ) S n+ 1 2 2
Ti .

This implies that the energy is positive provided that

1 2 C ∆t hi ≤ 1 c 1 - k c
. In other words, using quasiuniformity of the mesh, we conclude to the following CFL condition:

∆t h ≤ M c -k c b ,
with M a positive constant independent of the mesh and discretization parameters, for which the energy is positive.

Proposition 5. The modified energy verifies a discrete energy principle.

Proof. Using the scheme with the adequate test functions P n+1 + P n and S n+ 1 2 + S n-1 2 in respectively the first and second equation of ( 14), one obtains:

ε r P n+1 2 Ti -ε r P n 2 Ti + k c b × (S n+ 3 2 -S n+ 1 2 ), P n+1 + P n Ti = ∆t -S n+ 1 2 , ∇ × (P n+1 + P n ) Ti + ∆t S n+ 1 2 * × n, P n+1 + P n ∂Ti , µ r S n+ 3 2 2 Ti -µ r S n+ 1 2 2 Ti - k c b × (P n+1 -P n ), S n+ 3 2 + S n+ 1 2 Ti = ∆t P n+1 , ∇ × (S n+ 3 2 + S n+ 1 2 ) Ti -∆t P n+1 * × n, S n+ 3 2 + S n+ 1 2 ∂Ti . (16)
We add the two equations and rewrite them using the properties of the mixed product:

ε r P n+1 2 Ti + 2 k c b × S n+ 3 2 , P n+1 Ti + µ r S n+ 3 2 2 Ti -ε r P n 2 Ti -2 k c b × S n+ 1 2 , P n Ti -µ r S n+ 1 2 2 Ti = -∆t S n+ 1 2 , ∇ × (P n+1 + P n ) Ti + ∆t S n+ 1 2 * × n, P n+1 + P n ∂Ti + ∆t P n+1 , ∇ × (S n+ 3 2 + S n+ 1 2 ) Ti -∆t P n+1 * × n, S n+ 3 2 + S n+ 1 2 ∂Ti . ( 17
)
The right hand side of the last inequality can be decomposed in the sum of two elements I 1 i and I 2 i , with

I 1 i := -∆t S n+ 1 2 , ∇ × P n+1 Ti + ∆t S n+ 1 2 * × n, P n+1 ∂Ti + ∆t P n+1 , ∇ × S n+ 1 2 Ti -∆t P n+1 * × n, S n+ 1 2 ∂Ti , and
I 2 i := -∆t S n+ 1 2 , ∇ × P n Ti + ∆t S n+ 1 2 * × n, P n ∂Ti + ∆t P n+1 , ∇ × S n+ 3 2 Ti -∆t P n+1 * × n, S n+ 3 2 ∂Ti .
Using Green's formula, one has:

I 1 i := ∆t S n+ 1 2 , P n+1 × n ∂Ti + ∆t S n+ 1 2 *
× n, P n+1 ∂Ti -∆t P n+1 * × n, S n+ 1 2 ∂Ti , Thus:

I 1 i := -∆t 1 2 P n+1 + × n, S n+ 1 2 ∂Ti -∆t 1 2 P n+1 × n, S n+ 1 2 + ∂Ti ,

Finally,

i

I 1 i = 0.
For I 2 i , we wish to keep the symmetry of the expression and thus write:

I 2 i := - 1 2 ∆t S n+ 1 2 , ∇ × P n Ti - 1 2 ∆t ∇ × S n+ 1 2 , P n Ti + 1 2 ∆t S n+ 1 2 , P n × n ∂Ti + ∆t S n+ 1 2 * × n, P n ∂Ti + 1 2 ∆t P n+1 , ∇ × S n+ 3 2 Ti + 1 2 ∆t ∇ × P n+1 , S n+ 3 2 Ti - 1 2 ∆t P n+1 , S n+ 3 2 × n Ti -∆t P n+1 * × n, S n+ 3 2 ∂Ti . (18)
This gives:

I 2 i := - 1 2 ∆t S n+ 1 2 , ∇ × P n Ti + ∇ × S n+ 1 2 , P n Ti + ∆t S n+ 1 2 + × n, P n ∂Ti + 1 2 ∆t P n+1 , ∇ × S n+ 3 2 Ti + ∇ × P n+1 , S n+ 3 2 Ti -∆t P n+1 + × n, S n+ 3 2 ∂Ti . (19)
Summing over all the elements, taking into account all the contributions, we find:

ε r P n+1 2 + 2 k c b × S n+ 3 2 , P n+1 + µ r S n+ 3 2 2 -ε r P n 2 -2 k c b × S n+ 1 2 , P n -µ r S n+ 1 2 2 = i I 2 i (20) Thus E n+ 3 2 h -E n+ 1 2 h = 0. ( 21
)
The last results gives, in combination to the CFL condition of Proposition 4, the stability of the scheme.

The same kind of analysis can be adapted to take into account for the upwind fluxes (the latter have only negative contributions in the energy principle), and to more complicated schemes such as explicit Runge-Kutta type time discretization scheme. This could be obtained by adapting the arguments of this paper and those of [START_REF] Lanteri | Analysis of a generalized dispersive model coupled to a DGTD method with application to nanophotonics[END_REF]. However, this would end up in a lengthy, fastidious and technical proof that we choose not to reproduce here. In the simulations, we will mostly use explicit LSRK schemes with fully upwind fluxes (as mentioned earlier) and we will verify empirically the accuracy of the above CFL condition. Remark 6. As a side note, we would like to notice that one could be tempted to discretize system (10) as follows with a LSRK scheme:

∂ t P i ∂ t S i n+1 = M -1 i -K i × S i + l∈Vi S il (S * × n il ) K i × P i -l∈Vi S il (P * × n il ) n + -k × ∂ t P i k × ∂ t S i n ,
the additional time-derivative term being stored at each time-step and acting a source on the right-hand side for the computation of the next time-derivative of the fields. Although this method leads to a stable and valid discretization, the stability constraints with increasing angles are much worse than that obtained with time integration of the scheme (10) using LSRK4, leading to unbearably long computation times.

We now specify the type of sources and boundary conditions that we are using.

Incident field

Obviously, the field transform method is restricted to source fields respecting the PBC in the case of normal incidence. Here, we cover the well-known case of the wide-band plane wave source in a background material. It reads for respectively the electric and magnetic fields

E inc (r, t) = E 0 sin ω 0 t - k • r c b exp   - t -k•r c b 2 2σ 2    , H inc = k × E inc Z b , k • E 0 = 0, ( 22
)
where ω 0 is the central frequency of the pulse, σ controls the width of the Gaussian (and hence the spectral range), and Z b is the optical impedance of the background material. The third equation expresses that the polarization and the wavevector must be orthogonal. Applying transformation (3) to the frequency-domain version of system (22) and transforming back to time-domain is straightforward:

P inc (r, t) = P 0 sin ω 0 t - k⊥ • r c b exp   - t -k⊥ •r c b 2 2σ 2    , S inc = k × P inc Z b , k • P 0 = 0.
Care must be taken that the transformation does not affect the second and third relations, where the full wavevector still appears.

CFS-PML

We now turn to the treatment of the exterior boundary. As mentioned earlier, we have to artificially truncate the domain in the non-periodic directions. We could have used absorbing boundary conditions, but we rather choose to add a PML layer at the artificial boundary that will absorb the scattered field. Initially designed by Kuzuoglu et al [START_REF] Kuzuoglu | Frequency dependence of the constitutive parameters of causal perfectly matched layers[END_REF], the first adaptation of CFS-PMLs to the DGTD framework was recently proposed in [START_REF] König | Discontinuous Galerkin methods in nanophotonics[END_REF]. These PMLs rely on a complex stretching of the spatial coordinates:

∂ ∂k → 1 s k (ω) ∂ ∂k with s k (ω) = 1 - σ k iω -α k and k ∈ {x, y, z} . (23)
In the latter expression, σ k represents the damping rate of the PML, while α k is the actual frequency shift. Initially, a real stretch was also present, although it was removed here for the sake of simplicity. The reader is referred to [START_REF] König | Discontinuous Galerkin methods in nanophotonics[END_REF] for additional details about the CFS-PML. On the PML layer, the modified Maxwell's equations for the scattered field under this stretch are given by:

ε r ∂E ∂t = ∇ × H -G E + -G E -, ∂G E + ∂t = σ + ∇ + × H -(α + + σ + ) G E + , ∂G E - ∂t = -σ -∇ -× H -(α -+ σ -) G E -, µ r ∂H ∂t = -∇ × E -G H + -G H -, ∂G H + ∂t = -σ + ∇ + × E -(α + + σ + ) G H + , ∂G H - ∂t = σ -∇ -× E -(α -+ σ -) G H -, (24)
where the following notations are used:

∇ + × U =   ∂ y U z ∂ z U x ∂ x U y   , ∇ -× U =   ∂ z U y ∂ x U z ∂ y U x   , σ + =   σ y σ z σ x   , σ -=   σ z σ x σ y   ,
and similarly for α + and α -. The notation designates the element-wise multiplication for vectors, while G U + and G U -are additional fields required for the resolution of the field damping within the PML region. Applying the field transformation (3) to the full system of Maxwell's equations with CFS-PMLs would lead to a very large, tangled system with 18 equations, thus considerably enlarging the size of matrix Q. To get around this problem, it is however possible to apply the coordinate stretching (23) directly to the transformed system (4). Since the modifications brought by the stretching only impact the spatial derivatives, the standard derivation can be followed as in [START_REF] König | Discontinuous Galerkin methods in nanophotonics[END_REF], and the final continuous system of equations reads (again for the scattered field ):

ε r ∂P ∂t + k c b × ∂S ∂t = ∇ × S -G P + -G P -, ∂G P + ∂t = σ + ∇ + × S -(α + + σ + ) G P + , ∂G P - ∂t = -σ -∇ -× S -(α -+ σ -) G P -, µ r ∂S ∂t - k c b × ∂P ∂t = -∇ × P -G S + -G S -, ∂G S + ∂t = -σ + ∇ + × P -(α + + σ + ) G S + , ∂G S - ∂t = σ -∇ -× P -(α -+ σ -) G S -. (25)
After discretization, the CFS-PML equations can be solved as they are with the standard Maxwell's scheme (see [START_REF] Viquerat | Simulation of electromagnetic waves propagation in nano-optics with a high-order discontinuous Galerkin time-domain method[END_REF]), while the damping fields G U + and G U -(for U ∈ {P, S}) only need to be incorporated to the RHS of the P and S update equations. The reader is referred to [START_REF] König | Discontinuous Galerkin methods in nanophotonics[END_REF] and [START_REF] Viquerat | Simulation of electromagnetic waves propagation in nano-optics with a high-order discontinuous Galerkin time-domain method[END_REF] for the details of the DG discretization of CFS-PMLs.

Total Field/Scattered Field (TF/SF) formulation.

The combination of the absorbing PML layer and the incident field is dealt with by using a classical TF/SF strategy that exploits the possibility of splitting the contribution of the fields in their incident and scattered parts. In this manner, the computational domain is split into two parts: one containing the main scatterer where the total field is computed and the other complementary part where only the scattered field is computed (thus including the CFS PML zone). The incident field is in consequence imposed through the TF/SF interface that delimits these two zones (see e.g. [START_REF] Viquerat | Simulation of electromagnetic waves propagation in nano-optics with a high-order discontinuous Galerkin time-domain method[END_REF] for more details) and figure 1 for the whole picture.

Incorporation of dispersion models

At the frequency range considered in nanophotonics, the impact of dispersion can not be neglected in metallic media. The electrons don't instantaneously react to the applied electric field. This is rendered via frequency dependent permittivity laws, and thus frequency dependent speed of propagation of a given incident wave. In the temporal domain, this behavior is expressed via the introduction of one (or possibly multiple) polarization current(s). In the simplest version of a dispersion model (that is the framework used in this work), the evolution of the polarization currents is driven by auxiliary ODE's that are linearly coupled to Maxwell's equation as source currents.

For the sake of clarity, we recall the most classical version of such dispersion models called the generalized dispersion model.

µ 0 ∂ t H = -∇ × E, ε 0 ε s ∂ t E = ∇ × H -J f -J b , ∂ t J f = -γJ f + ω 2 P ε 0 E, J b = J 0 + i∈L1 J i + i∈L2 J i , J 0 = σ + i∈L2 J i E, J i = a i E -b i P i , ∀i ∈ L 1 , ∂ t P i = J i , ∀i ∈ L 1 , ∂ t J i = (c i -d i f i )E -f i J i -e i P i , ∀i ∈ L 2 , ∂ t P i = d i E + J i , ∀i ∈ L 2 . (26)
The fields J f and J b are the polarization currents of respectively the free and bound electrons. The parameters γ, ω 2 P and σ are given physical parameters. The bound polarization decouples into several polarization currents. Their evolution is given using the set of parameters (a i , b i ) i∈L1 , (c i , d i , e i , f i ) i∈L2 , that are given parameters, used to design the model according to experimental data. The classical Drude model, as referred to in this paper, corresponds to the case where J b = 0.

We refer the reader to [START_REF] Viquerat | Simulation of electromagnetic waves propagation in nano-optics with a high-order discontinuous Galerkin time-domain method[END_REF][START_REF] Lanteri | Analysis of a generalized dispersive model coupled to a DGTD method with application to nanophotonics[END_REF] for more details.

In the sequel for realistic numerical simulations, we will rely on such models. The extension of the field-transform technique is straightforward in this case.

Numerical results

We now concentrate on the numerical results. We will first provide some useful observables. Then, we illustrate the behavior inferred from the theoretical results and then tackle several nanophotonic test cases.

Observables

A few observables are usually at the heart of the analysis of periodic structures. In this section, we shortly review some of them, and see how they transform in the case of oblique incidence.

Reflection and tranmission

The most common quantities of interest for periodic structures are certainly the reflectance R and the transmittance T , which are frequency-dependent quantities that provide informations about the capability of a structure to reflect or let power flow through it. They are usually computed from the time-averaged Poynting vector:

π = 1 2 Ê × Ĥ * , (27)
where we denote by • the Fourier transform. Integrated over proper monitoring surfaces, this leads to the knowledge of R and T :

R (ω) = Si π sca • n Si π inc • n , T (ω) = So π tot • n Si π inc • n
, where the subscripts refer to scattered, incident and total field.

Thanks to the conjugate in expression (27), applying the field transform (3) simply leads to:

π = 1 2 P × Ŝ * . (28)

Diffraction efficiency

A particularly important class of periodic structures is that of gratings, which have the ability to diffract incident light in a discrete set of beams, each one propagating in a specific direction. The functioning of a grating is usually characterized by the computation of its diffraction efficiencies η n,m , which describe how the transmitted (or reflected) power divides between the existing diffraction orders [START_REF] Palmer | Diffraction Grating Handbook[END_REF]. Assuming a doubly periodic domain with lateral sizes a x and a y , the propagative scattered far field at a fixed height z 0 above the grating can be decomposed on a Fourier basis:

Êsca (r, ω) = (n,m)∈Z 2 e n,m (ω) exp i k sca x x + k sca y y = (n,m)∈Z 2 e n,m (ω) exp [i ((k x + k n ) x + (k y + k m ) y)] , (29)
where we defined k n = 2πn ax , k m = 2πm ay . In expression (29), each couple (n, m) contributes to the sum under the condition that:

k sca z = |k| 2 -(k x + k n ) 2 -(k y + k m ) 2 ∈ R, (30)
meaning that only propagative modes can contribute to the far field. The complex Fourier coefficients are defined by:

e n,m (ω) = 1 a x a y ax ay Ê(r, ω)| z0 exp [-i ((k x + k n ) x + (k y + k m ) y)] dxdy. ( 31
)
Analogous steps apply for the H field, leading to a similar expression for h n,m (ω). For each propagative mode, an associated Poynting vector can be computed:

π n,m (ω) = 1 2 e n,m (ω) × h * n,m (ω) , (32)
from which the mode-wise diffraction efficiency can be defined:

η n,m (ω) = a x a y π n,m (ω).

(33)

The diffraction efficiency is related to the reflectance as follows, under condition (30):

R(ω) = (n,m)∈Z 2 η n,m (ω), (34)
As can be seen from expressions (35) and (33), the definition of η n,m remains unchanged under transformation (3):

π n,m (ω) = 1 2 p n,m (ω) × s * n,m (ω) , (35)
where p n,m and s n,m are defined transparently from expression (31).

Time-step reduction

As a first validation of the results of section 2.2, we consider a bi-periodic domain of lateral size 75 nm filled with vacuum, terminated with CFS-PMLs regions at the top and bottom. A broadband plane wave with significant intensity in the range λ ∈ [150, 600] nm is injected through a Total-Field/Scattered-Field (TF/SF) interface, and the quantity ∆t = ∆t(θ=0 • ) ∆t(θ) for θ ∈ [0 • , 90 • [ is computed. The scheme used is LSRK4 with upwind fluxes (α = 1). The results, obtained for piecewise P 3 polynomial approximation, are shown in figure 5 along with λ x + . Although there is not a perfect match, the general trend confirms our previous analysis. The deviation observed for large incidence angles can be attributed to the dissipation induced by the fully upwind scheme that stabilizes the numerical algorithm.

Critical angle with dielectric background

We now consider a slab of vacuum (c r = 1) embedded in a dielectric background for which c b = 1 √ 2 , i.e. ε b = 2. In this case, the expected critical angle is θ c = 45 • . The rest of the configuration is similar to that of last section. In figure 6, we plot the figure of merit ∆t : as before, although the match is not perfect, the trend of the numerical results decently follows the theory. As it was suggested in section 2.2, no stable time-step can be found for θ > θ c , confirming the previous analysis.

Bi-periodic gold slab

We consider the case of a bi-periodic gold slab of thickness 50 nm embedded in vacuum and illuminated by a wide-band plane wave pulse with TM polarization. As before, the lateral size of the domain is 75 nm, and the incident field is injected through a TF/SF interface. The gold is described by a simple Drude model, and the wavelength range of interest is [150, 600] nm. In this problem, the observable of interest is the reflectance, which analytical solution can be computed with a transfer matrix method. Regarding the DGTD discretization, a P 3 polynomial approximation is used along with a fourth-order scheme in time with upwind fluxes (α = 0), leading to an overall fourth-order accuracy in both space and time. Results are shown in figure 7. A perfect match is observed between the numerical and the exact solution, both for normal and oblique incidences.

λ (nm) R(λ) 0°2 0°5 0°F
igure 7: Reflectance of a bi-periodic 50 nm thick gold slab under normal and oblique incidence computed with third-order polynomial basis (full line) compared to exact solution (dots).

Dielectric color filter

We now consider optical color filters composed of silicon nanodisks deposited on a glass substrate and embedded in PMMA, as presented in [START_REF] Park | Structural color filters enabled by a dielectric metasurface incorporating hydrogenated amorphous silicon nanodisks[END_REF] (see figure 8). In this contribution, the authors present a set of substractive filters, for which they compute and measure the transmission in normal incidence. In this section, we reproduce the results for a magenta filter under normal incidence, before exploring its response under oblique incidence using the field transform technique, for various angles and polarizations. Results are shown in figure 9. As can be seen, at oblique incidences, the filter response is altered: in TM polarization, the filtered wavelength is blue-shifted for θ = 20 • , while the quality factor of the filter is seriously decreased once 40 • of incidence is reached. In TE polarization, the filtered wavelength is red-shifted for θ = 20 • with new resonances appearing for longer wavelengths. For θ = 40 • , the filter efficiency is seriously compromised. As an illustration, we show in figure 10 a field map of the Fourier transform of E x at the operating wavelength of the color filter, illuminated by a monochromatic plane wave with incidence angle θ = 30 • . For clarity, the fields are not represented in the PML regions (at the bottom and the top of the domain). For the magenta filter, the height of the silicon cylinder is 80 nm, and its diameter is 130 nm. The lateral size of the domain is 330 nm. At the bottom, the glass is embedded directly in the PML region, thus mimicking an infinitely thick glass layer.

Two dimensional silicon grating

We now consider a two dimensional silicon grating on a silicon substrate of infinite depth as depicted in figure 11. The lateral periodic unit cell length is d G = 700 nm, the height of the grating h G = 350 nm, and the substrate infinitely extended in the negative z direction. For silicon, we used the material parameters from [START_REF] Green | Self-consistent optical parameters of intrinsic silicon at 300 K including temperature coefficients[END_REF]. As for the previous examples, the incident field is injected through a TF/SF interface, on which the reflectance R(λ) and diffraction efficiencies η n,m (λ) are computed. Results of the computation for normal (θ = 0 • ) and oblique (θ = 20 • ) incidences are shown in figures 13 and 12. In normal incidence, it is visible that the normal mode η 0,0 is the only contributor to the reflectance above 500 nm. At lower wavelengths, a set of four symmetric high-order diffraction modes appear. As could be expected, the situation is different at oblique incidence, where high-order diffraction modes appear at larger wavelengths. Only two symmetric modes appear below 600 nm, while a higher-order asymmetric mode is present below 480 nm. For both polarizations, the filtering function is altered, either by diminishing its efficiency, or by shifting the filtered wavelength.

Conclusions

In this paper, we complemented the work of [START_REF] Miller | A discontinuous Galerkin time domain framework for periodic structures subject to oblique excitation[END_REF] on the use of the field transformation technique to handle oblique incidence with regular PBC on 3D periodic structures in the DGTD framework. The transformed equations were analyzed with the method of characteristics and energy techniques, highlighting several limitations that were then confirmed through numerical experiments. The DGTD discretization of the transformed system was recalled, and the stability of the semi-discrete formulation was proven. Then, insights were given about the effects of the field transformation on the incident field, the observables, and the of CFS-PMLs in this context. Our implementation was validated on a textbook case for which exact solution exists, before being applied to the study of two real-life cases: (i) the robustness of a siliconbased color filter various incidence angles, and (ii) the computation of the efficiency of high-order diffraction modes in a silicon grating. ). The reflectance coefficient is shown in gray dots, while the relevant diffraction efficiencies are shown in blue, cyan and orange. As could be expected, the tilted incidence angle gives rise to different high-order diffraction angles than for the normal incidence case (figure 12). Only two symmetric modes appear below 600 nm, and an additional asymmetric mode is also present below 480 nm.

Figure 1 :

 1 Figure1: 2D configurations of normal and oblique incidence with PBC. The domain is periodic in the lateral directions. and is infinite in the vertical direction. Here the notation PMLs at top and bottom already refers to the artificial numerical truncation of the domain that will be used in the sequel. The incident plane wave is oriented toward the bottom of the domain. Numerically, this will be imposed at an artificial interface called the Total-Field/Scattered-Field interface (see later for details). In normal incidence, the PBC applies naturally, and at each discretization point matching information is available from the periodic neighbor point (green dots). In oblique incidence, a causality issue appears, since the matching information for the left green dot is located in a zone that has not been illuminated by the incident field yet.

Figure 2 :

 2 Figure 2: Definition of (θ, φ) angles in the cartesian coordinate system.

Figure 3 :

 3 Figure 3: Normalized characteristic speeds of system (5) for φ = 0 and c = 1.

Figure 4 :

 4 Figure 4: Plot of the critical value θc as a function of c in the range [0, 1[. For c > 1, θ can take all values in [0, 90 • [.

Figure 5 :

 5 Figure 5: Relative maximal stable timestep ∆t = • ) ∆t(θ)for θ ∈ [0 • , 90 • [ compared to λ x + in vacuum.

Figure 6 :

 6 Figure 6: Relative maximal stable timestep ∆t = ∆t(θ=0 • ) ∆t(θ)for θ ∈ [0 • , 45 • [ compared to λ x + for a slab of vacuum embedded between two dielectric half-spaces.

Figure 8 :

 8 Figure8: Silicon-based color filter embedded in PMMA, deposited on a glass substrate. For the magenta filter, the height of the silicon cylinder is 80 nm, and its diameter is 130 nm. The lateral size of the domain is 330 nm. At the bottom, the glass is embedded directly in the PML region, thus mimicking an infinitely thick glass layer.

Figure 9 :

 9 Figure9: Transmission of the magenta color filter at normal and oblique incidences, for both TE and TM polarizations. For both polarizations, the filtering function is altered, either by diminishing its efficiency, or by shifting the filtered wavelength.

Figure 10 :

 10 Figure 10: Field map of the Fourier transform of Ex for the color filter setup. The device is illuminated from below by a monochromaticplane wave with incidence angle θ = 30 • .

Figure 11 :

 11 Figure 11: Computational mesh of the grating unit cell. Gray: periodic boundary faces; pink: TF/SF interface; red: silicon grating on infinite substrate; the remaining domain is air. Both the silicon substrate and the air are truncated by PMLs (not shown here) in the ẑ direction.

Figure 12 : 1 Figure 13 :

 12113 Figure12: Diffraction efficiencies of the grating at normal incidence (θ = 0 • ). The reflectance coefficient is shown in gray dots, while the relevant diffraction efficiencies are shown in dark and light blue. As can be seen, the normal mode η 0,0 is the only contributor to the reflectance above 500 nm, where a set of four symmetric modes appear.


whether they are published or not. The documents may come   L'archive ouverte pluridisciplinaire

Related research

► Unfamiliarity -L2 listeners typically find it difficult to understand speakers with an unfamiliar, often non-standard accent (Major et al. 2005, Ockey et

Results

Percentage of correctly transcribed syllables in a recording of a speaker with the accent being studied.

Statistical significance (p=0.0199) between Pretest and after both film and sentence training. (No statistical significance for each step individually.)

  al. 2016, inter alia) ► Intelligibility vs. comprehensibility -Speech can be rated highly for intelligibility but low for comprehensibility.(Derwing & Munro 1997) -Intelligibility -Whether an utterance can be understood -Comprehensibility -The amount of perceptual and cognitive effort required to understand an utterance ► Accent intelligibility training -Intelligibility can be improved through training for L1 listeners -ESL teachers acquire this skill naturally (Gass & Veronis 1984) -L1 listeners can be trained in a lab setting (Floccia et al. 2009 inter alia) ► Kasparek (2008) -Sentence-level training showed better results for L1 listeners of an unfamiliar accent than word-level training. Helping L2 Learners Improve their Understanding of Different Accents An in-class experiment adapting methods used in laboratory settings Experiment: Adaptation to an L2(+) classroom setting ► Context --Half-semester C1 class on accents in English (16 hours). -Mix of L1 languages (French, Arabic, Chinese, Spanish, German). Included in study: eight men, three women. ► Two-phase experiment --1) Watching a film mostly in Hindi-speaking Indian accent (Amu) -2) Two 30-minute transcription sessions with sentences borrowed from film Discussion and Future Work ► Working with full sentences seems to help improve intelligibility for L2(+) students in classroom setting ► Challenges: Maintaining student interest; similar level of difficulty in tests ► Future work: Testing sentence-level training before watching film -Removing influence of film in testing for significance -"Rewarding " students with film