Ludovic Marigo-Lombart

Christophe Viallon

Alexandre Rumeau

Alexandre Arnoult

Aurelie Lecestre

Laurent Mazenq

Ayad Ghannam

Hugo Thienpont

Krassimir Panajotov

Guilhem Almuneau
email: almuneau@laas.fr.

Aurélie Lecestre

A Lecestre

Electro-Absorption Modulator vertically integrated on a VCSEL: microstrip-based high-speed electrical injection on top of a BCB layer

Keywords: High-speed modulation, Electro-absorption modulator, VCSEL, BCB

teaching and research institutions in France or abroad, or from public or private research centers.

I. INTRODUCTION

I N order to increase the modulation bandwidth, the vertical integration of an electro-optic modulator onto a VCSEL (EOM-VCSEL) has already been demonstrated to be a competitive approach [START_REF] Barron | Millimeter-wave asymmetric fabryperot modulators[END_REF]- [START_REF] Yakimov | Concept of feedback-free highfrequency loss modulation in detuned duo-cavity vertical cavity surfaceemitting laser[END_REF]. In that case, a double mesa structure is required with three contacts: one on the top for applying a voltage to the EOM, one middle contact for a shared ground, and one at the bottom for the VCSEL CW current injection.

This approach requires : 1) planarization with a low permittivity layer to decrease the parasitic capacitances; 2) characterization of dielectric properties of such a layer up to frequencies as high as 100 GHz; 3) development of a MS access line to avoid the losses induced by the electromagnetic (EM) coupling of the signal with the doped substrate. In this paper, we address all the above issues in the case of an EAM-VCSEL, in which a voltage-driven electro-absorber placed above modulates the output light of the steady-biased VCSEL.

Indeed, all the groups that have already worked on vertical integration of an EOM onto a VCSEL used highly doped substrate for the current injection in the VCSEL without considering, to our knowledge, the induced EM losses for the modulator. In [START_REF] Paraskevopoulos | Ultra-high-bandwidth (> 35GHz) electrooptically-modulated VCSEL[END_REF], an optical modulation bandwidth exceeding 35 GHz was demonstrated using p-doped substrate, intracavity injection scheme for the modulator and CPW access line. Chen et al [START_REF] Chen | Multilevel Amplitude Modulation Using A Composite-Resonator Vertical-Cavity Laser[END_REF] worked with a Ground-Signal-Signal-Ground (GSSG) access line that presents the advantage of providing all the access on the same level. However, the substrate doping remains an issue and the fabrication process complexity is considerably increased. Finally, Yakimov et al. [START_REF] Yakimov | Concept of feedback-free highfrequency loss modulation in detuned duo-cavity vertical cavity surfaceemitting laser[END_REF] used a doped substrate with a CPW access line. Although a p-doped substrate has been chosen for the EAM-VCSEL presented in this paper, we show that an undoped substrate could potentially improve the performance of a CPW-based device at the cost of a more complex process for VCSEL fabrication that will not be discussed here.

In order to increase the EOM electrical frequency bandwidth, a low-permittivity material is required for separating the modulator electrode pads. To this aim, metal air bridge can be used but this approach suffers from a complicated process flow. A thick dielectric layer can also be implemented but it is difficult to shape it by photolithography and etching while ensuring precise pattern alignment. Finally, polyimides are very good candidates and are most often used to overcome these last issues. Compared to SiO x or Si x N y , polyimides present lower dielectric constant [START_REF] Ganrou | Stress-buffer and passivation processes for Si and GaAs IC's and passive components using photosensitive BCB: Process technology and reliability data[END_REF]. We develop in this paper a novel approach to planarize Bis-benzocyclobutene (BCB) on top of semiconductor devices by using a mechanical press benefiting from a nanoimprint tool to apply a high pressure and temperature and spread the BCB and improve its flatness. BCB is nowadays of high interest in many applications such as wafer bonding, 3D interconnect processes or high-speed devices in MEMS, silicon and photonics applications. Indeed, it presents good mechanical properties after hard cure, low moisture and low permittivity ( r ). Furthermore, BCB is more easily definable by spin-coating and photolithography compared to dielectric layers that require plasma or chemical etching, which can significantly damage the semiconductor topmost layers. These BCB properties, despite its relatively low thermal conductivity, make it very attractive for microwave optoelectronics. Indeed, BCB is extensively used in all high speed semiconductor based optical devices such as Distributed Feedback Lasers, Vertical-Cavity Surface-Emitting Lasers (VCSEL) and electro-optical modulators.

The precise knowledge of BCB permittivity and losses for frequencies up to 100 GHz is crucial for predicting and optimizing the influence of the RF electrical access line on device performance. The dielectric constant of BCB provided by the manufacturer is 2.65 at 1 MHz and 2.5 in the 1-10 GHz frequency range, and within a temperature range from 30 to 200 • C, whatever the BCB solution [8]. The few charaterization results available from the literature report r values in the range of 2.5-2.7, depending on the measurement technique and the used BCB solution [START_REF] Dahlgren | An integrated millimeterwave BCB patch antenna HEMT receiver[END_REF]- [START_REF] Al-Omari | Dielectric characteristics of spincoated dielectric films using on-wafer parallel-plate capacitors at microwave frequencies[END_REF]. Furthermore, BCB has been mostly characterized at frequencies below 40 GHz. In [START_REF] Sakai | A millimeter-wave flip-chip IC using microbump bonding technology[END_REF], a value of 2.7 has been extracted up to 60 GHz by fitting the Sparameters of a 50 Ω MS line involving a 26 µm-thick BCB layer. In [START_REF] Costanzo | Benzocyclobutene as Substrate Material for Planar Millimeter-Wave Structures: Dielectric Characterization and Application[END_REF], 2.65 was measured in the range 11-65 GHz from a CPW. In this paper, the BCB permittivity and the losses are characterized up to 100 GHz.

In order to provide carrier injection in VCSELs, a CPW access line is most commonly used [START_REF] Larsson | High-speed VCSELs for short reach communication[END_REF], [START_REF] Moser | Energy Efficiency of Directly Modulated Oxide-Confined High Bit Rate 850-nm VCSELs for Optical Interconnects[END_REF]. To avoid EM coupling with the highly doped bottom n-contact layer close to the active region, the n-layer is etched away under the pad line which constitutes an additional process step. However, this approach is not possible for EAM-VCSEL, in which a highly doped substrate must be used as a third electrode. Therefore, we develop in this paper an optimized MS access line. We study the influence of the substrate doping and compare the high frequency properties of the MS line configuration to that of CPW access line on N, P and Not Intentionally Doped (NID) substrates with several test structures.

II. BCB INNOVATIVE PLANARIZATION PROCESS A. BCB technological implementation

As described in [START_REF] Burdeaux | Benzocyclobutene (BCB) dielectrics for the fabrication of high density, thin film multichip modules[END_REF], BCB can be simply spin-coated with as good degree of planarization (DOP) as 95 %. This optimum DOP is achieved for device pitch of about 400-500 µm but it worsens when the distance between components decreases down to 150-200 µm, which is commonly used for VCSEL manufacturing. The double mesa structure of the EAM-VCSEL is also a major hurdle to get the aimed surface flatness. Furthermore, the BCB process requires an additional photolithography step with a tight alignment in order to open an access down to the semiconductor surface.

To overcome all these issues, it is possible to planarize the BCB at the height level of the mesa structure. To do so, chemical mechanical polishing (CMP) can be used to lower the BCB thickness after thermal curing but this process is difficult and may lead to unsatisfying results [START_REF] Modafe | Embedded benzocyclobutene in silicon: An integrated fabrication process for electrical and thermal isolation in MEMS[END_REF] and possible contamination. Alternatively, one can etch the BCB in a controlled way by plasma reactive etching. However, in the case of non-perfect BCB planarization, a slight swelling is created after spin-coating above the mesa edge and is replicated after dry etching. This can cause a short-circuit as shown in the topmost mesa corner.

Finally, a mechanical pressing techniques using nanoimprint tool can be advantageously and reliably used for polymer planarization as previously demonstrated in [START_REF] Shih | CON-TACT R planarization process of spin-on dielectrics for device fabrication[END_REF] and [START_REF] Abada | Uniform fabrication of thick SU-8 patterns on small-sized wafers for micro-optics applications[END_REF]. After obtaining a sufficient BCB flatness, Yang et al. [START_REF]Reactive ion etching characteristics of partially-cured benzocyclobutene[END_REF] reported a study on reactive ion etching (RIE) conditions for BCB etching with SF 6 -based gas mixtures. A stop-etch layer (for example, Si 3 N 4 ) can be used to accurately match the aimed etch depth [START_REF] Demir | Self-Aligning Planarization and Passivation for Integration Applications in III-V Semiconductor Devices[END_REF].

B. Self-aligned BCB planarization process

We present in this section an innovative process combining a step of mechanical pressing in order to flatten the BCB surface with a step of reproducible dry etching of the BCB excess film, while keeping the surface smooth enough for high-frequency low-loss injection. Prior to the BCB coating, a SiO x layer is deposited by inductively coupled plasma chemical vapour deposition (ICPCVD) to improve BCB adhesion [START_REF] Marigo-Lombart | Self-aligned BCB planarization method for high-frequency signal injection in a VCSEL with an integrated modulator[END_REF]. After a spin-coating with ramping-up speed and a soft bake of BCB 4026-46 resist, a layer thickness of 10.6 µm is obtained. After the soft-bake, done at 80 • C, it is still possible to reflow the BCB layer by heating it to temperature of 100 • C for 6 minutes. We also demonstrate that neither higher temperature nor longer heating are necessary. Taking advantage of this reflow, a mechanical press is used, with a nano-imprint equipment under 1.72 MPa and 100 • C for 6 minutes in order to perfectly flatten the BCB layer.

After the planarization and before the dry etch, it is mandatory to anneal the BCB by hardcure treatment at 250 • C for one hour. In such a way most of the contained solvent is removed and the BCB mechanical strength and electrical properties are improved.

The final step of the BCB planarization is the thickness readjustment in order to strictly outcrop the mesa heights, thus resulting in a perfect planarized configuration whatever the spacing between etched mesas. This etching step is applied on the overall surface of the sample without any masking materials. Different conditions have been tested for the dry etching of BCB layer with various gases, pressure and, ICP and RF powers. By inductively coupled plasma reactive ion etching (ICP-RIE), the lowest roughness with controllable etching rate is achieved with a O 2 /SF 6 ratio of 98/2 at 80 mTorr and ICP and bias powers of 500 W and 40 W, respectively. More details on this process can be found in [START_REF] Marigo-Lombart | Self-aligned BCB planarization method for high-frequency signal injection in a VCSEL with an integrated modulator[END_REF].

The BCB surface roughness directly influences the quality of the Au top contact and so the high-frequency performance. A slight roughness is however required to ensure a good adhesion of the metal on the BCB, and is naturally created by the dry etch process. We monitor in-situ the reflected 670 nm laser signal: its average amplitude providing the etch depth, while the amplitude oscillations providing the etched surface roughness. AFM measurements show a roughness of only 10 nm on the metal layer deposited on partially etched BCB.

III. BCB MICROWAVE CHARACTERIZATION

Most commonly, the extraction of the dielectric microwave characteristics of insulating materials is done by inserting the studied material into a waveguide or a transverse electric and magnetic (TEM) transmission line. The complex dielectric constant and permeability are then retrieved from the measured scattering matrix as described in [START_REF] Weir | Automatic measurement of complex dielectric constant and permeability at microwave frequencies[END_REF].

In our case, only a thin layer of BCB can be deposited onto a substrate of different nature, which makes it virtually impossible to design a purely TEM structure. We therefore chose to focus on the quasi-TEM MS transmission line for its ease of fabrication: the BCB layer is placed between the signal path and the lower ground plane. Since the EM wave propagates partly in the air and partly in the material to be characterized, measurements give only access to an effective permittivity. The proven analytical models that are available for the case of the MS line are then used to retrieve the relative dielectric permittivity from the measured effective one.

A. Wafer fabrication

Only a few steps are required to fabricate MS line structures for BCB characterization. First, a ground plane located under the BCB is deposited on the whole sample. Then, a BCB layer is spin-coated and an additional photolithography step is employed to form openings used to electrically connect the ground plane to the ground pads located on the top of the BCB layer. Finally, after a heating treatment, a lift-off is performed to fabricate the RF pads, the vias and the transmission line.

In the case of an EAM vertically integrated onto a VCSEL, the first mesa is 3 µm thick and another etch step is required to reach the active region of the VCSEL and to provide lateral electrical and optical confinement. That corresponds to a global mesa thickness of 6 µm.

Several MS lines with lengths of 200, 300, 700, 1200 µm and a 14 µm-width signal track have been drawn, as done in [START_REF] Ghannam | Dielectric microwave characterization of the SU-8 thick resin used in an above IC process[END_REF]. These lines present a characteristic impedance of 53 Ω using Cyclotene TM manufacturer standard permittivity of 2.65. Coplanar squared RF pads with 60 µm widths are added at the MS lines terminals to enable Ground-Signal-Ground (GSG) RF probe contacts.

B. Access pads de-embedding

The first step prior to the extraction of BCB dielectric characteristics is to remove the RF pads contribution from raw measured data. Usually, an electrical equivalent lumped model of the pad is retrieved from dedicated Short, Open and/or Thru dummies structures [START_REF] Song | A de-embedding technique for interconnects[END_REF]. Open and Short patterns are used to calculate the equivalent parallel capacitance and series inductance of the pad, respectively. But these patterns are becoming problematic with increasing frequency due to the non-ideality of the corresponding terminations [START_REF] Ito | A simple through-only de-embedding method for on-wafer S-parameter measurements up to 110 GHz[END_REF]. Such problems are avoided if the Open and Short devices are replaced by only one thru pattern but this technique introduces new issues [START_REF] Ito | A simple through-only de-embedding method for on-wafer S-parameter measurements up to 110 GHz[END_REF], [START_REF] Tretiakov | A new on-wafer de-embedding technique for on-chip rf transmission line interconnect characterization[END_REF]. Indeed, since the pad is modeled by lumped elements, all these methods are valid only for access dimensions much smaller than the wavelength (typ. λ g / 20) of the highest test frequency. A short thru pattern is hence desirable but the close proximity of the signal pads and the probes introduces a cross-talk that leads to an incorrect pad modeling [START_REF] Tretiakov | A new on-wafer de-embedding technique for on-chip rf transmission line interconnect characterization[END_REF]. The technique proposed here avoids all these issues. Two transmission lines of different lengths are required instead of Open/Short or Thru dummies structures. The method described below borrows an important part of [START_REF] Saavedra-Gómez | On-wafer CMOS transistors de-embedding method using two transmission lines of different lengths[END_REF] but differs from it by the assumptions used to simplify the problem to make it consistent and by the subsequent nature of the model chosen for the pad. In [START_REF] Saavedra-Gómez | On-wafer CMOS transistors de-embedding method using two transmission lines of different lengths[END_REF], the access pads are modeled by a distributed reciprocal and symmetric (S 11 = S 22 ) two-port network. Here, the symmetric condition is impaired because of the large step in the widths between the pad to the MS line that leads to S 11 = S 22 , especially at high frequencies. Removing this symmetry leads to an under-determined problem that is handled by re-introducing a simplified lumped pad model similar to the one used in the Open/Short or Thru de-embedding methods.

T A T B T l1 T A T B T l2 (a) T A T B Z Z Y Y (b)
Both transmission lines display the same characteristic impedance Z c and propagation constant γ but differ by the lengths, written as l 1 and l 2 , respectively. The surrounding signal pads are modeled by the lumped networks drawn in Fig. 1b. Using these parameters, T A and T B , T l1 and T l2 are expressed as follows

T A = 1 Z Y Y Z + 1 , T B = Y Z + 1 Z Y 1 (2)
T l1 = cosh(γl 1 ) Z c sinh(γl 1 ) 1
Zc sinh(γl 1 ) cosh(γl 1 )

T l2 = cosh(γl 2 ) Z c sinh(γl 2 ) 1 Zc sinh(γl 2 ) cosh(γl 2 ) (3)
From Fig. 1a, we write

M 1 = T A T l1 T B (4)
M 2 = T A T l2 T B . (5)
Expressions ( 4) and ( 5) are re-writen as follows

M 1 T B -1 = T A T l1 (6)
M 2 T B -1 = T A T l2 . (7)
Since transmission lines are symmetric devices, m 11 = m 22 and p 11 = p 22 are observed experimentally. Using these relations and expanding ( 6) and [START_REF] Ganrou | Stress-buffer and passivation processes for Si and GaAs IC's and passive components using photosensitive BCB: Process technology and reliability data[END_REF], the following expressions are derived

m 12 = [cosh(γl 1 ) + m 11 ] Z 1 + Y Z + Z c 1 + Y Z sinh(γl 1 ) (8) m 21 = [cosh(γl 1 ) + m 11 ]Y + 1 + Y Z Z c sinh(γl 1 ) (9)
p 12 = [cosh(γl 2 ) + p 11 ] Z 1 + Y Z + Z c 1 + Y Z sinh(γl 2 ) (10) p 21 = [cosh(γl 2 ) + p 11 ]Y + 1 + Y Z Z c sinh(γl 2 ). (11)
After defining the matrix A as

A = cosh(γl 1 ) + m 11 sinh(γl 1 ) cosh(γl 2 ) + p 11 sinh(γl 2 ) , (12)
the expressions (8) to [START_REF] Janezic | Permittivity characterization of low-k thin films from transmission-line measurements[END_REF] are rewritten and combined in a matrix form using [START_REF] Costanzo | Benzocyclobutene as Substrate Material for Planar Millimeter-Wave Structures: Dielectric Characterization and Application[END_REF] to get the final equation

Y Z 1+Y Z 1+Y Z Zc Zc 1+Y Z = A -1 • m 21 m 12 p 21 p 12 . (13)
Equation ( 13) can be solved for Y , Z and Z c knowing M 1 and M 2 and the related propagation constant γ. γ can be extracted beforehand from M 1 and M 2 , i.e. by using the method proposed in [START_REF] Reynoso-Hernández | Computation of the wave propagation constant gamma in broadband uniform millimeter wave transmission line[END_REF].

C. Extraction of the electrical characteristics of BCB

First, the complex effective permittivity eff is calculated from the de-embedded measurement data of one of the two transmission lines. The transmission-reflection method introduced in [START_REF] Weir | Automatic measurement of complex dielectric constant and permeability at microwave frequencies[END_REF] and modified in [START_REF] Ghannam | Dielectric microwave characterization of the SU-8 thick resin used in an above IC process[END_REF] is used to extract eff . The relative permittivity r , as well as the loss tangent term tan δ are then determined either by a direct computation carried out by using the analytic equations derived from Bahl's model (method 1) [START_REF] Bahl | Simple and accurate formulas for a microstrip with finite strip thickness[END_REF], or by using an iterative algorithm based on the Jensen-Hammerstad model (method 2) [START_REF] Hammerstad | Accurate Models for Microstrip Computer-Aided Design[END_REF]. The choice between these two methods is dictated by the validity range of each model regarding the MS line geometry. It mainly depends on the thickness over height t/h ratio of the analyzed MS line. These two techniques, as well as the form factor validity range have been investigated in [START_REF] Ghannam | Dielectric microwave characterization of the SU-8 thick resin used in an above IC process[END_REF]. Accurately measured geometric data of the line are also required to insure the best possible accuracy of the extraction of the BCB electrical characteristics.

The measured S 21 parameters of all MS lines are shown in Fig. 2a. A vector network analyser Anritsu Vectorstar, its Broadband Test Set 3738A and a probe station were employed to perform measurements up to 110 GHz. 150 µm pitch Cascade Infinity GSG probes were also used.

The BCB permittivity r and tan δ extracted by using the two previously described methods reach the same outcome, as seen in Fig. 2b. At low frequency, a permittivity of 2.6 is measured. It decreases to 2.5 at 70 GHz. Past 70 GHz, we go beyond the validity range of the lumped element approximation used to model the pads, which probably induces the change of slope observed in the r vs frequency curve. tan δ is defined as the ratio between imaginary and real parts of the complex permittivity ( / ). Its values are comprised between 0.01 and 0.03 in the 0.05-70 GHz range.

The method developed here to extract the BCB permittivity from MS lines measurements will be reused in Section V-B to calculate the access line contribution on the measured impedance of the modulator.

IV. HIGH-FREQUENCY INJECTION DESIGN

The measured BCB dielectric data are now used to design the pads and the access line and to minimize the HF injection losses. In this section, two designs are compared: the CPW commonly used by the VCSEL community, and the MS access line proposed here.

A. CPW access line 1) Pad design: we based our first approach on a standard coplanar structure where both the ground plane and the access line are located on the top of the BCB layer as shown in Fig. 3a. The ground plane surrounds the signal pad to promote ground equipotentiality and avoid slot-line propagation mode. For clarity, only heights are on scale, lateral dimensions have been shortened. Here, the signal pad is tapered to decrease the reflection coefficient at its extremity and constitutes a simple line to the modulator top contact. The distances shown in Fig. 3a and Fig. 4a between the ground plane and the signal pad and line are named g 1 and g 2 , respectively. It is necessary to bring the ground plane very close to the access line in order to confine the EM wave in the CPW structure and limit its interaction with the underlying substrate. With a BCB thickness of 6 µm, this confinement is verified using an EM simulator (ANSYS-HFSS) by setting g 1 and g 2 to 15 µm and 5 µm, respectively. These last dimensions allow to decrease as much as possible the penetration depth of the magnetic field in the doped substrate while considering the fabrication constraints for the metal lift-off.

Furthermore, the BCB thickness is decreased to a few micrometers around the EAM mesa in order to reach the top VCSEL mesa and thus the intermediate contact. This last dimension is crucial and presents a very difficult process step since it is necessary to control the sidewall inclination to reach the VCSEL contact. Indeed, the VCSEL mesa diameter should not exceed 60-70 µm to avoid problems related to long oxidation times in order to define the lateral confinement.

2) Process: Three different substrate-dopings have been used: n-doped, p-doped and Not Intentionally Doped (NID). Above the deposited BCB layer prepared as described in II-B, e-beam evaporated Ti/Au 50 / 200 nm bi-layer is used to form the pads.

3) Impact of the substrate doping with the CPW design: CPW thru structures formed by interconnecting the inner end of two EAM access lines, as shown in Fig. 4a, have been fabricated on the three substrate types and the aforementioned BCB layer. Measurements have been carried out using the same setup as in III-C. In order to understand the influence of substrate doping on the CPW access losses, the equivalent electrical Π-model shown in Fig. 4b is extracted from the measured S ij data for NID, n and p-doped substrates. The values of R s , L s , R p and C p are plotted in Fig. 5. The CPW on the NID substrates exhibits a very small R s resistance around 1 Ω below 60 GHz, 2 to 4 times lower than with doped substrates. With doped substrates, an important part of the return current goes through the substrate. The induced negative magnetic coupling increases the resistance R s and slightly lowers the inductance L s . A highly doped substrate also acts as a poor metal layer that shields the CPW line and increases its equivalent capacitor C p compared to NID substrate. Undoped substrates are preferred in recent high-speed VCSELs for all these reasons [START_REF] Larsson | High-speed VCSELs for short reach communication[END_REF]. Finally, the performance of CPW access lines is strongly dependent on the substrate conductivity. When a heavily doped substrate is used, the induced access resistance R s and capacitance C p are increased and may limit the maximum operating frequency of the modulator. This result motivated the study of a MS access line.

B. MS access line

The MS line can avoid the previously described substrateinduced effects if a ground plane supporting the return current is placed above the substrate. The electric contact required between the EAM and the VCSEL in order to apply bias voltage is used as a MS ground plane. In such way, the RF performance of this MS configuration becomes independent on the substrate doping. Here, we compare it with the previously analyzed CPW access line. The shared ground contact is now moved under the 6 µm-thick BCB layer, while the access line remains on top of it. In such a way, the complicated angled photolithography step is no longer required, which substantially facilitates the fabrication process of the double mesa structure. As it is still necessary to have GSG contacts on the top of the chip to connect the RF probe, so the BCB layer is opened to bring the contacts to the ground pads as seen in Fig. 3b. The dimensions of this MS line design will be discussed in the next section. We now compare in Fig. 6 the CPW and the MS access lines with the same lengths (130 µm). Despite the losses brought by the substrate conductivity in the case of CPW access, the p-doped subtrate is used for this comparison since it minimizes the optical losses between the EAM and the VCSEL devices. The CPW and the MS line widths are 19 and 14 µm, respectively. The metal thickness is 0.3 µm for the CPW. Thicknesses of 0.3 µm / 0.3 µm and 0.6 µm / 0.7 µm are studied for the signal/ground layers of the MS line.

In Fig. 6, the EM simulation data are in good agreement with our measurements on CPW. The CPW exhibits higher access resistance despite the larger line width because of the magnetically-induced substrate losses as previously explained. Doubling the MS line thickness has a very limited effect especially at the highest frequencies because of skin depth limitations. In both cases, the access inductance is comparable but the capacitance is much lower with a MS line thanks to its narrower signal line. Finally, even though the gain with the MS design is not so significant, this geometry is preferred due to an easier fabrication process.

V. HIGH-FREQUENCY OPERATION OF THE EAM

A. Device description

The GaAs-based modulator studied in this paper is composed of an undoped active region formed by multiple quantum wells (MQWs) playing the role of electro-absorber. The MQWs are embedded between two distributed Bragg reflectors (DBR) to provide an asymmetric Fabry-Perot structure. More details about the EAM structure can be found in [START_REF] Marigo-Lombart | Vertical electro-absorption modulator design and its integration in a VCSEL[END_REF]. No compositional grading is added at the DBR interfaces for the sake of simplicity of the epitaxial growth.

The modulation is achieved by applying an electric field perpendicular to the structure through quantum confined Stark effect (QCSE). QCSE modifies the band structure of the MQWs and thus the envelope function of the electrons and holes. In the case of the monolithic integration of the EAM onto a VCSEL, the VCSEL is CW powered by injecting current between the bottom electrical contact and the middle ground contact. The modulation frequency (RF) signal is applied between the top electrode of the EAM and the middle ground contact as shown in Fig. 3. The MS access line and the open vias to transfer the ground contact on top are presented in Fig. 3, too.

EAM samples are processed in the same way as for complete EAM-VCSEL devices, meaning that we fabricate the same double mesa structure to keep the same BCB thickness on the standalone EAM and on the EAM-VCSEL.

The operation of the standalone EAM device with a MS access line has already been demonstrated in [START_REF] Marigo-Lombart | High frequency characterization of a vertical electro-absorption modulator for data communications[END_REF], using a RF synthesiser source and a spectrum analyzer to measure the photodiode response. As shown in Fig. 7, modulation bandwidth above 28 GHz has been achieved, limited by the photodiode responsivity.

B. Equivalent circuit

In order to further understand the EAM behavior at highfrequency, its electrical response is extracted from the measured reflection coefficient. The EAM is biased through 50 to 150 µm MS lines terminated by GSG pads. To de-embed the signal pad and the access line contributions, two lines displaying the same geometrical characteristics but with different lengths are processed on the same wafer. The electrical characteristics of the signal pad and the line are recovered by using the method described in section III-B. Let T be the chain matrix comprizing the matrices T A of the pad and T l of the access transmission line

T = T A T l = t 11 t 12 t 21 t 22 . (14)
The de-embedded modulator impedance Z mod is deduced from the measured impedance Z meas by using the expression Z mod = t 22 Z meas -t 12 t 11 -t 21 Z meas [START_REF] Moser | Energy Efficiency of Directly Modulated Oxide-Confined High Bit Rate 850-nm VCSELs for Optical Interconnects[END_REF] We are now able to get a better insight in the access element contributions to the overall impedance of the EAM. The EAM acts like a capacitance at low frequency, because of the 485 nm-thick active region. Moreover, its impedance is more important for smaller diameter as seen in Fig. 8. At high frequency, the impedances of the structures with two different diameters tend to get closer to each other, which reflects a modification of the electrical behavior of the EAM. It can now be considered as a resistance in series with a capacitance more than a single capacitance. We can notice after the de-embedding, that the impedance of the EAM alone is higher than the one of the global circuit. Indeed, the presence of the pad and the line between the probe extremity and the modulator decreases the efficiency of the electrical injection in the device. The difference between the two curves, with and without de-embedding, is furthermore increased when the diameter is smaller (higher impedance) since the pad and line contributions, compared to the EAM impedance, are more important. This clearly emphasizes the relevance of the de-embedding process to know precisely the electrical behavior of the EAM.

VI. CONCLUSION

In conclusion, thanks to an original and reliable planarization technological process, we have extracted a BCB permittivity of 2.56 up to 70 GHz with a variation of ±2.5 % in this frequency range. Based on these characterizations, we designed two RF access geometries in order to optimize the high frequency modulation voltage actually applied to the EAM. Since working with doped substrate is mandatory for monolithically integrated EAM-VCSEL, CPW appears unsuitable due to the substrate-induced losses. It is shown that MS access line configuration makes easier the fabrication of a double mesa while being providing superior RF performances compared to CPW. Thanks to the method we developed, we have de-embedded the pad and the access line and extracted the impedance of the standalone EAM. This impedance could be drastically improved with DBR graded interfaces, in the same way as in a current-injected standard VCSEL. Finally, reducing further the EAM mesa diameter would considerably increase the electrical bandwidth of the EAM-VCSEL device.

Fig. 1 . 1 = m 11 m 12 m 21 m 22 and M 2 = p 11 p 12 p

 11212 Fig. 1. (a) Transmission lines L 1 and L 2 with related chain matrices. (b) Lumped models of pads A and B used to calculate T A and T B . The two lines L 1 and L 2 and the terminating pads are shown in Fig. 1a. If the ABCD formalism is used, the measured scattering parameters of lines L 1 and L 2 are converted into chain parameters, marked M 1 and M 2 , respectively M 1 = m 11 m 12 m 21 m 22 and M 2 = p 11 p 12 p 21 p 22 .(1)

FrequencyFig. 2 .

 2 Fig. 2. (a) Measured transmission parameter of MS lines of 200 (black), 300 (red), 700 (blue) and 1200 µm (magenta). (b) Extracted BCB data.

Fig. 3 .

 3 Fig. 3. (a) CPW and (b) MS designs for EAM-VCSEL device.

Fig. 4 . 4 Fig. 5 .

 445 Fig. 4. (a) design of CPW thru device and (b) Π equivalent lumped network used to model this device.

Fig. 6 .

 6 Fig.6. Rs, Ls, Rp and Cs elements of measured (red) and simulated (black) CPW access compared to simulations of MS access with thin metal deposition (blue) and thicker metal (dark blue).

Fig. 7 .

 7 Fig. 7. Frequency modulation response of EAM devices of different diameters : 18 µm (black), 23 µm (blue) and 26 µm (red).

Fig. 8 .

 8 Fig. 8. Raw ( ) and deembedded ( ) EAM impedance Z mod for 19 µm (red) and 25 µm (black) devices.

The authors gratefully acknowledge the technological support of RENAT-ECH (French Network of Technology Platforms) within LAAS-CNRS cleanroom infrastructure and the financial support by the Methusalem foundation, Belgium.

Manuscript received xxxxx; revised xxxx.


Basile De Loynes
email: basile.deloynes@ensai.fr

Fabien Navarro
email: fabien.navarro@ensai.fr

Baptiste Olivier
email: baptiste.olivier@orange.fr

‡ Baptiste

Data-driven Thresholding in Denoising with Spectral Graph Wavelet Transform

Keywords: Spectral Graph Theory, Denoising, Stein Unbiased Risk Estimation, Spectral Graph Wavelet Transform, Tight Frame, Variance Estimation

This paper is devoted to adaptive signal denoising in the context of Graph Signal Processing (GSP) using Spectral Graph Wavelet Transform (SGWT). This issue is addressed via a data-driven thresholding process in the transformed domain by optimizing the parameters in the sense of the Mean Square Error (MSE) using the Stein's Unbiased Risk Estimator (SURE). The SGWT considered is built upon a partition of unity making the transform semiorthogonal so that the optimization can be performed in the transformed domain. However, since the SGWT is over-complete, the divergence term in the SURE needs to be computed in the context of correlated noise. Two thresholding strategies called coordinatewise and block thresholding process are investigated. For each of them, the SURE is derived for a whole family of elementary thresholding functions among which the soft threshold and the James-Stein threshold. This multi-scales analysis shows better performance than the most recent methods from the literature. That is illustrated numerically for a series of signals on different graphs.

Introduction

The emerging field of Graph Signal Processing (GSP) aims to bridge the gap between signal processing and spectral graph theory (see for instance [START_REF] Chung | Spectral graph theory[END_REF]; [START_REF] Belkin | Towards a theoretical foundation for Laplacian-based manifold methods[END_REF] and references therein). One objective is to generalize fundamental analysis operations from regular grid signals to irregular structures as graphs. There is an extensive literature on GSP, in particular we refer the reader to [START_REF] Shuman | The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains[END_REF] for an introduction to this field and [START_REF] Ortega | Graph signal processing: Overview, challenges, and applications[END_REF] for an overview of recent developments, challenges and applications. As a matter of fact, GSP have already been applied in machine/deep learning: convolutional neural networks (CNN) on graphs [START_REF] Bruna | Spectral networks and locally connected networks on graphs[END_REF]; [START_REF] Henaff | Deep convolutional networks on graph-structured data[END_REF]; [START_REF] Defferrard | Convolutional neural networks on graphs with fast localized spectral filtering[END_REF], semi-supervised classification with graph CNN [START_REF] Kipf | Semi-supervised classification with graph convolutional networks[END_REF]; [START_REF] Hamilton | Inductive representation learning on large graphs[END_REF], community detection [START_REF] Tremblay | Graph wavelets for multiscale community mining[END_REF], to name just a few. In the context of GSP, the authors of [START_REF] Coifman | Diffusion maps[END_REF]; [START_REF] Gavish | Multiscale wavelets on trees, graphs and high dimensional data: Theory and applications to semi supervised learning[END_REF]; [START_REF] Hammond | Wavelets on graphs via spectral graph theory[END_REF] have developed wavelet transforms on graphs. More specifically, in [START_REF] Hammond | Wavelets on graphs via spectral graph theory[END_REF] a fairly general construction of a frame enjoying the usual properties of standard wavelets is developed: each vector of the frame is localized both in the graph domain and the spectral domain. The transform associated with this frame is named Spectral Graph Wavelet Transform (SGWT).

Many studies based on SGWT (or some variants) explore the denoising performance of this approach using different strategies [START_REF] Leonardi | Tight wavelet frames on multislice graphs[END_REF]; [START_REF] Onuki | Graph signal denoising via trilateral filter on graph spectral domain[END_REF]; [START_REF] Wang | Trend filtering on graphs[END_REF]; [START_REF] Deutsch | Manifold denoising based on spectral graph wavelets[END_REF]; [START_REF] Irion | Efficient approximation and denoising of graph signals using the multiscale basis dictionaries[END_REF]; [START_REF] Dong | Multiscale representation of surfaces by tight wavelet frames with applications to denoising[END_REF]; [START_REF] Göbel | Construction of tight frames on graphs and application to denoising[END_REF] from signal adapted tight frames to regularization method.

The denoising approach chosen in this paper involves several thresholding processes in the transformed domain of the wavelet coefficients. Actually, this can be seen as an extension to SGWT of the methodology of [START_REF] Donoho | Adapting to unknown smoothness via wavelet shrinkage[END_REF]; [START_REF] Cai | Adaptive wavelet estimation: a block thresholding and oracle inequality approach[END_REF]. With this approach, the main challenge is the efficient calibration of the parameters minimizing the MSE risk in a complete data-driven way. Recently, in the setting of discrete wavelets transform on a regular grid-the so-called regular case-the Stein's unbiased risk estimate (SURE) has proven to be a powerful tool for signal/image restoration [START_REF] Luisier | A new SURE approach to image denoising: interscale orthonormal wavelet thresholding[END_REF]; [START_REF] Pesquet | A SURE approach for digital signal/image deconvolution problems[END_REF] ;[START_REF] Vaiter | Local behavior of sparse analysis regularization: applications to risk estimation[END_REF]. Based on the Stein's lemma, this estimator acts as a proxy for the MSE which cannot be computed in practice since the original signal is unknown. In this paper, the SURE is explicitly computed for an arbitrary thresholding process in Theorem 1 and for correlated noise in the graph domain in Corollary 1. Also, let us point out that contrary to the regular wavelet transform, the SGWT is no longer orthogonal so that a white Gaussian noise in the graph domain is transformed in a correlated noise. Consequently, the divergence term of the resulting SURE involves the covariance of the transformed noise making the numerical evaluation less simple than in the regular case. Afterward, the SURE is specified to the case of coordinatewise and block thresholding. The latter is inspired by image denoising problems for which a Stein risk estimator has been proposed in [START_REF] Peyré | Adaptive structured block sparsity via dyadic partitioning[END_REF] to tune both the block-sparsity structure and the threshold. A similar selection strategy has been developed by [START_REF] Navarro | Adaptive parameter selection for block wavelet-thresholding deconvolution[END_REF] in the context of deconvolution. The R package gasper which implements the method introduced in this paper is available on github1 (de [START_REF] De Loynes | Gasper: Graph signal processing in r[END_REF] as well as the scripts to reproduce the results presented2 . Many denoising methods consist in shrinking the coefficients of a sparse representation of the considered signals in a certain transformed domain. In [START_REF] Talebi | Global image denoising[END_REF], the aim is to attenuate the coefficients in the Fourier domain corresponding to deletion of the high frequencies contained in the noisy signal. The selection of the parameters is done by minimizing the SURE estimating the MSE in a similar way to our methodology. Due to the limited ability of ideal lowpass graph filter to separate the low-frequency noise, an adaptive weighted graph filter is proposed in [START_REF] Chen | Image denoising with adaptive weighted graph filtering[END_REF]. These two graphs filtering methods can be approximately mimicked with the SWGT based one by removing the finest scales and an appropriate choice of the thresholding. However, the localization properties of SGWT both in space and frequency allow to better process signals with heterogeneous regularity in the initial domain. Intuitively, the thresholding process in the SGWT transformed domain eliminates the high frequency part of the noise where the signal is locally regular and tends to keep the peaks of the original signal whereas the Fourier transform tends to smooth them. This comes at the cost of a redundant frame in which the transformed noise is no longer stationary neither decorrelated.

Regarding the regularization method implemented in [START_REF] Onuki | Graph signal denoising via trilateral filter on graph spectral domain[END_REF], the regularization parameter is also selected optimizing an MSE proxy based on a similar argument. Nonetheless, beyond the fact that the philosophy is different (regularization versus thresholding), one stress that the empirical risk bias is explicitly determined while the MSE estimation in [START_REF] Onuki | Graph signal denoising via trilateral filter on graph spectral domain[END_REF] is only validated numerically. Another penalization method is given in [START_REF] Wang | Trend filtering on graphs[END_REF] that extends the approach from [START_REF] Tibshirani | The solution path of the generalized lasso[END_REF] within the framework of graphs. For this method, the divergence term is computed explicitly; this gives rise to a datadriven parameter selection method so that this approach is an interesting concurrent to our methodology.

The paper is organized as follows. Section 2 introduces the notation and briefly reviews the notions of tight frame and SGWT of [START_REF] Hammond | Wavelets on graphs via spectral graph theory[END_REF]. Section 3 is devoted to denoising and the SURE estimator for generic thresholding process in the transformed domain. Then, the SURE is specified in the cases of coordinatewise, block thresholding processes and for correlated noise in the graph domain. In Section 4 numerical comparisons with the classical Wiener filter (oracle version) and the trend filtering introduced in [START_REF] Wang | Trend filtering on graphs[END_REF] for denoising are discussed. Several signals and graphs, including examples from real datasets, are considered. For these experimental results, the construction of the frame follows [START_REF] Göbel | Construction of tight frames on graphs and application to denoising[END_REF]. In terms of denoising performance, other tight frames such as spectrum adapted and/or signal adapted tight frames from [START_REF] Shuman | Spectrum-adapted tight graph wavelet and vertex-frequency frames[END_REF] and [START_REF] Behjat | Signal-adapted tight frames on graphs[END_REF] might give better results. Still, Theorem 1 actually applies to any tight frame and the question of exhibiting the most efficient one is beyond the scope of the paper (see [START_REF] Shuman | Localized spectral graph filter frames: A unifying framework, survey of design considerations, and numerical comparison[END_REF] for a comprehensive survey).

2 Spectral Graph Wavelet Transform

Graphs, Frames and Tight Frames

Let G be an undirected weighted graph, with set of vertices V , and weights (w ij ) i,j∈V satisfying w ij = w ji for i, j ∈ V . The size of the graph is the number of nodes n = |V |. The (unnormalized) graph Laplacian matrix L ∈ R V ×V associated with G is the symmetric matrix defined as L = D -W , where W is the matrix of weights with coefficients (w ij ) i,j∈V , and D the diagonal matrix

with diagonal coefficients D ii = j∈V w ij . A signal f on the graph G is a function f : V → R. Let F = {r i } i∈I be a frame of vectors of R V , that is a family of vectors in R V such that there exist A, B > 0 satisfying for all f ∈ R V A f 2 2 ≤ i∈I | f, r i | 2 ≤ B f 2 2 . (1)
The linear map

T F : R V → R I defined for f ∈ R V by T F f = ( f, r i ) i∈I is called the analysis operator.
The synthesis operator is the adjoint of T F : namely, it is the linear map T * F : R I → R V defined for a vector of coefficients (c i ) i∈I by T * F (c i ) i∈I = i∈I c i r i . As a frame is in particular a generating family of R V , a signal f ∈ R V can be recovered from its coefficients T F f with the help of the synthesis operator.

Construction of Tight Frames

A frame F is said to be tight if A = B = 1 in Equation (1)-the latter is then termed the Parseval identity. From now on, the frames considered are supposed to be tight. Let us recall the generic construction of such a frame (c.f. [START_REF] Kereta | Monte carlo wavelets: a randomized approach to frame discretization[END_REF] for instance).

Since L is self-adjoint, it admits the spectral decomposition L = λ χ , • χ , where λ 1 ≥ λ 2 ≥ • • • ≥ λ n = 0 denote the (ordered) eigenvalues of the matrix L, and (χ ) 1≤ ≤n are the associated normalized and pairwise orthogonal eigenvectors. Then, for any function ρ : sp(L) → R defined on the spectrum sp(L) of matrix L, the functional calculus formula reads ρ(L) = ρ(λ ) χ , • χ . A finite collection (ψ j ) j=0,...,J is a finite partition of unity on the compact

[0, λ 1 ] if ψ j : [0, λ 1 ] → [0, 1] for all j ∈ J and ∀λ ∈ [0, λ 1 ],
J j=0 ψ j (λ) = 1. Given a finite partition of unity (ψ j ) j=0,...,J , the Parseval identity implies that the following set of vectors is a tight frame:

F = ψ j (L)δ i , j = 0, . . . , J, i ∈ V .
Also, following [START_REF] Leonardi | Tight wavelet frames on multislice graphs[END_REF]; [START_REF] Göbel | Construction of tight frames on graphs and application to denoising[END_REF], a partition of unity can be easily defined as follows: let ω : R + → [0, 1] be some function with support in [0, 1], satisfying

ω ≡ 1 on [0, b -1 ],
for some b > 1, and set

ψ 0 (x) = ω(x) and ψ j (x) = ω(b -j x) -ω(b -j+1 x) for j = 1, . . . , J, where J = log λ 1 log b + 2.

Discrete SGWT Associated with a Partition of Unity

Let (ψ j ) j=0,...,J be a partition of unity of [0,

λ 1 ]. The SGWT of a signal f ∈ R V is given by Wf = ψ 0 (L)f T , . . . , ψ J (L)f T T ∈ R n(J+1) .
The adjoint linear transformation W * of W is:

W * η T 0 , η T 1 , . . . , η T J T = j≥0 ψ j (L)η j .
The tightness of the underlying frame implies that

W * W = Id R V so that a signal f ∈ R V
can be recovered by applying W * to its wavelet coefficients ((Wf ) i ) i=1,...,n(J+1) ∈ R n(J+1) (see [START_REF] Hammond | Wavelets on graphs via spectral graph theory[END_REF]).

Adaptive Denoising with SGWT

Let f ∈ R V be some signal on a graph G and ξ be an n-dimensional Gaussian vector distributed as N (0, σ 2 Id). The aim of denoising is to recover the unknown signal f from the observed noisy version f = f + ξ. Basically, our denoising procedure will consist of three steps: (1) compute the

SGWT transform W f ∈ R n(J+1) ;
(2) apply a given thresholding operator h(•) to the coefficients W f ; (3) apply the inverse SGWT transform to obtain an estimation f of the original signal.

Here, the main challenge in denoising consists in choosing a suitable thresholding operator with respect to the noisy signal f and the underlying graph. The performance measure in the sequel will be the MSE between the original signal f and the denoised signal f : f -f 2 2 . First, it is worth noting that the Parseval identity allows direct optimization in the transformed domain of wavelet coefficients. Secondly, in practice, obviously the original signal remains unknown.

To overcome this difficulty, the MSE is generally substituted with the Stein's Unbiased Risk Estimator which no longer depends on the original signals (see [START_REF] Donoho | Adapting to unknown smoothness via wavelet shrinkage[END_REF] for instance). Nonetheless, contrary to the usual wavelet transform, the white noise ξ is mapped onto a correlated Gaussian noise. In the next section, the SURE is derived taking into account these correlations.

The SURE Estimator in the Transformed Domain

By linearity, the denoising problem f = f + ξ is transferred to the denoising problem [START_REF] Göbel | Construction of tight frames on graphs and application to denoising[END_REF] for instance). Consequently, its spectral decomposition reads WW * = U ΣU * with U a unitary matrix of R n(J+1) and Σ = Id R n 0 0 0 . A thresholding process is a map h : R n(J+1) → R n(J+1) . Typically, the map h is a coordinatewise or a block shrinkage in applications. The following result extending the SURE's expression to correlated noise is based on the Stein's lemma in [START_REF] Stein | Estimation of the mean of a multivariate normal distribution[END_REF] in which h is assumed to be weakly differentiable. One refers the reader to [START_REF] Stein | Estimation of the mean of a multivariate normal distribution[END_REF] for the precise definition.

F = F + Ξ with Ξ ∼ N (0, σ 2 WW * ), F = W f and F = Wf . It is worth noting that WW * is an orthogonal projector of rank n in R n(J+1) (see
Theorem 1 (h-SURE). Let h be a weakly differentiable thresholding process for the denoising problem F = F + Ξ. Then the theoretical MSE is given by

E h( F ) -F 2 = E   -nσ 2 + h( F ) -F 2 + 2 n(J+1) i,j=1 Cov(Ξ i , Ξ j )∂ j h i ( F )   , where h i is the i-th component of h.
It is worth noting that Cov(Ξ i , Ξ j ) = σ 2 (WW * ) i,j so that, as soon as the thresholding process h is specified and the noise variance σ 2 estimated, the SURE of h defined below can be completely computed from the noisy observations as in the regular case:

SURE(h) = -nσ 2 + h( F ) -F 2 + 2 n(J+1) i,j=1 Cov(Ξ i , Ξ j )∂ j h i ( F ).
Proof. The theoretical MSE can be rewritten as follows

E h( F ) -F 2 = E h( F ) -F 2 + E Ξ 2 + 2E h( F ) -F , Ξ .
The second term is equal to nσ 2 since almost surely Ξ 2 = U * Ξ 2 = P K U * Ξ 2 where K = ker(WW * ) ⊥ and P K the orthogonal projection onto K. Finally, setting g

(x) = h(x) -x, x ∈ R n(J+1) , it remains to compute the last term E g( F ), Ξ = E g(F + Ξ), Ξ where F is deterministic and Ξ ∼ N (0, σ 2 WW * ). A simple computation gives E g(F + Ξ), Ξ = n(J+1) i=1 E[g i (F + Ξ)Ξ i ] = n(J+1) i=1 Cov(g i (F + Ξ), Ξ i ).
Then, following [START_REF] Liu | Siegel's formula via Stein's identities[END_REF], each term in the sum above is given by

Cov(g i (F + Ξ), Ξ i ) = -nσ 2 + n(J+1) j=1 Cov(Ξ i , Ξ j )E[(∂ j h i )(F + Ξ)],
since Tr(σ 2 WW * ) = nσ 2 . This ends the proof of Theorem 1.

Coordinatewise Thresholding Process

For a coordinatewise thresholding process, the map h is of the form h(x) = (τ (x i , t i )) i=1,...,n(J+1) where (t i ) i=1,...,n(J+1) are the thresholds. In practice, we may choose τ (x, t) = x max{1t β |x| -β , 0} with β ≥ 1. The most popular choices are the soft thresholding (β = 1), the James-Stein thresholding (β = 2) and the hard thresholding (β = ∞). The latter will not be considered here since it does not lead to a sufficiently regular thresholding process for Theorem 1 to be applied.

For any β ∈ [1, ∞), the derivative ∂ j h i vanishes whereas

∂ i h i (F + Ξ) = 1 [t i ,∞) (| F i |) 1 + (β -1) t β | F i | β .
Consequently, the SURE associated with h is given by

SURE(h) = -nσ 2 + n(J+1) i=1 F 2 i 1 ∧ t β i | F i | β 2 + 2 n(J+1) i=1 V(Ξ i )1 [t i ,∞) (| F i |) 1 + (β -1)t β i | F i | β . (2)
The usual expression of the SURE is recovered from the identity above remarking that

V[Ξ i ]
are identically equal to σ 2 when the transformed noise is uncorrelated.

Let us notice that the coordinatewise soft thresholding (β = 1) satisfies an oracle inequality as shown in [START_REF] Göbel | Construction of tight frames on graphs and application to denoising[END_REF]. Similarly to the regular case, it states that up to a log factor, the soft thresholding estimator can mimic an oracle projection.

Optimization: Donoho and Johnstone's Trick

The SURE can be optimized in the same way as in the standard case using the Donoho and Johnstone's trick of [START_REF] Donoho | Adapting to unknown smoothness via wavelet shrinkage[END_REF] whose the justification is recalled below.

For the sake of simplicity, we first consider the case of the coordinatewise thresholding process with a uniform threshold: t i = t for all i = 1, . . . , n(J + 1). Denote by a 1 , . . . , a n(J+1)

the absolute values of the noisy wavelet coefficients | F i | in the increasing order. The trick comes from the observation that, on each interval (a k , a k+1 ), the last term of Equation ( 2) is nondecreasing whereas the second term is an increasing function of t. Consequently, the SURE hits its minimum at some value a k * , k * = 1, . . . , n(J + 1).

If the thresholds t i are no longer uniform but merely tied inside blocks with values t 1 , . . . , t L , the same trick is still valid: group the terms in the sums along the different parameters t 1 , . . . , t L and optimize each partial sum with respect to t k , k = 1, . . . , L.

Block Thresholding Process

In order to take advantage of the localization properties of SGWT and the regularity of the original signal, we may introduce block thesholding processes similar to [START_REF] Cai | Adaptive wavelet estimation: a block thresholding and oracle inequality approach[END_REF].

Consider a partition (B ) ∈L of {1, . . . , n(J + 1)} and set x 2 B = i∈B (x i ) 2 . In this case, the thesholding process h = (h i ) i=1,...,n(J+1) reads J+1) , and ∈ L : i ∈ B .

h i (x) = x i max 1 - t β x β B , 0 , x ∈ R n(
If i, j are in different blocks, then ∂ j h i vanishes. Additionally, if i, j are in B but i = j then

∂ j h i ( F ) = F i 1 [t ,∞) ( F B )βt β F j F -β-2 B , whereas
∂ i h i ( F ) = 1 [t ,∞) ( F B ) 1 -t β F -β B + βt β F 2 i F -β-2 B .
Consequently, a straightforward computation leads to

SURE(h) = -nσ 2 + ∈L 1 ∧ t β F β B 2 F 2 B + 2 ∈L 1 [t ,∞) ( F B )   1 - t β F β B i∈B V(Ξ i ) + βt β F β+2 B i,j∈B Cov(Ξ i , Ξ j ) F i F j  
Once again, for uncorrelated transformed noise, the usual expression easily follows from the identity above. Note also that the optimization of the SURE in this case requires more sophisticated techniques as the divergence term is no longer monotone.

Correlated Noise in the Graph Domain

The SURE can also be stated in the context of correlated noise at the cost of some prior information on the covariance structure. More precisely, in the denoising problem f = f + ξ with correlated noise, it is supposed that ξ ∼ N (0, Γ) for some covariance matrix Γ. The denoising problem reads in the transformed problem as F = F + Ξ with Ξ ∼ N (0, WΓW * ).

Corollary 1. Under the assumption of Theorem 1, the theoretical MSE is given by

E[ h( F ) -F 2 ] = E -Tr(WΓW * ) + h( F ) -F 2 + 2 n(J+1) i,j=1 Cov(Ξ i , Ξ j )∂ j h i ( F ) ,
where h i is the i-th component of h.

Let us point out that the parameters selection can be made without computing explicitly Tr(WΓW * ) since it does not depend on h-even though, the MSE estimate is obviously shifted by this quantity. Besides, the correlation structure Γ is actually hidden in the quantities Cov(Ξ i , Ξ j ), namely, for 1 ≤ i, j ≤ n(J + 1): Cov(Ξ i , Ξ j ) = WΓW * i,j . Consequently, computationally speaking, there is no additional burden compared to the white noise case.

Proof. The proof follows the lines of Theorem 1 with E[Ξ

2 i ] = WΓW * i,i , 1 ≤ i ≤ n(J + 1).
In applications, it is usually reasonable to assume some structure on the covariance matrix Γ reflecting the topology of the underlying graph. Typically, the noise on two given vertices may be correlated if those vertices are close enough in the graph. For example, let ξ 0 ∼ N (0, σ 2 Id),

we set ξ = ξ 0 + αW ξ 0 where W is the graph matrix of weights and α > 0 some tuning parameter describing the global intensity of the correlation. Then, it follows,

Γ = σ 2 (Id + 2αW + α 2 W W * ). (3)
Other choices are obviously possible.

Complexity

Regarding the space complexity, we need to store the frame and the weights appearing in the SURE for a cost of O(n 2 (J +1) 2 ). With given Laplacian eigendecomposition, the time complexity of the optimization of the SURE is (in average) of order O(n(J + 1) log(n(J + 1))) for the coordinate-wise estimator following [8]. For the block estimator, the use of a grid search is a limitation.

Numerical Results

This

The Minnesota Roads Graph

The Minnesota roads graph is a planar graph consisting of 2642 vertices and 6606 edges. Each vertex is described by its (x, y)-coordinates. The function ω chosen in the experiments is a piecewise linear function with support in [0, 1] and constant equal to 1 on [0, b -1 ] with b = 2.

From λ 1 ≈ 6.89, we deduce that the number of scales is J + 1 = 5-see Section 2. We compare the performance in terms of SNR (computed on the functions after reconstruction) for different denoising strategies, different noise levels and each synthetic signal f η,k . For each noise level σ = 0.005, 0.01, 0.02 and σ = 0.001, 0.002, 0.004, a sample of N = 10 white Gaussian noise is simulated and a global (G) versus level-dependent (LD) coordinatewise thresholding are performed with respect to the soft (β = 1) and James-Stein (β = 2) thresholding rules. For each strategy, we compare the average behavior of the SNR for parameter selected with the oracle (MSE β=1,2 obtained by minimizing the MSE using the original signal f ) and the SURE with known σ (SURE β=1,2 σ ). Also, the standard deviation on the sample is provided.

These results are first compared to the classical Wiener filter. More precisely, the Wiener filter consists of attenuating the Fourier coefficient F( f ) of f . Below, we only consider the oracle

linear attenuation F( f )[i]F(f )[i]/(F(f )[i] 2 + σ 2 )
. While this estimator is unrealistic since it depends on f , any Wiener filter has worse performance than this oracle. Table 1 is completed by the performance of the Wiener filter on each signal. Also, the theoretical value r inf of the oracle risk given in [START_REF] Mallat | A wavelet tour of signal processing[END_REF] is recalled for comparison purpose.

Our methodology is also compared to the so-called graph trend filtering (i.e. for k = 0, 1, 2) introduced in [START_REF] Wang | Trend filtering on graphs[END_REF]. The graph trend filtering is a regularization method with a penalty term involving the graph difference operator at a given order (see [START_REF] Wang | Trend filtering on graphs[END_REF]).

In the experiments, we make use of the matlab toolbox gtf3 provided by the authors of [START_REF] Wang | Trend filtering on graphs[END_REF]. 1 that our method performs better than the trend filtering motivating the use of multiscale analysis. This idea is confirmed by the comparison with the Wiener filter, in particular in the lower SNR regime that is for higher noise levels. Also, similarly to the regular case, numerical experiments shows that the James-Stein threshold (β = 2) is slightly more efficient than the soft threshold in particular for the global thresholding process.

Let us point out that there is no fundamental difference in terms of performance between the global and level dependent thresholding in this experiment. In fact, the level dependent thresholding always performs at least as good as the global one. Since the additional computational cost is acceptable, the level dependent thresholding appears to be a good choice without any further a priori knowledge.

The Facebook Graph

Here we examined and compare the denoising performance of level dependent SGWT thresholding (LD) against the trend filtering and Laplacian smoothing [START_REF] Smola | Kernels and regularization on graphs[END_REF] on a nonplanar graph considered in [START_REF] Wang | Trend filtering on graphs[END_REF]: the Facebook graph from the Stanford Network Analysis Project 4 . This undirected graph, collected from survey participants using this Facebook app, is composed of 4039 nodes representing Facebook users, and 88,234 edges representing friendships (see [START_REF] Leskovec | Learning to discover social circles in ego networks[END_REF] for more details). For signal f , we consider the different regularities used in [START_REF] Wang | Trend filtering on graphs[END_REF] as well as the same noise levels, for 5 realizations (see (Wang et al., 2016, Section 5.1) for more details). More precisely, we simply run the scripts provided in the matlab toolbox provided by the authors of [START_REF] Wang | Trend filtering on graphs[END_REF].

The results are shown in the Figure 2 (to be compared with (Wang et al., 2016, Figure 5 comparable results, LD globally provides better performance than trend filtering (whatever the value of k), especially at high noise levels where the maximal gain in terms of SNR is greater than 15dB and exceeds 5dB and 10dB respectively for the 5 highest noise levels for inhomo- The calibration of certain parameters has not been studied. However, the latter can considerably influence the performance of the SGWT. For the homogeneous random walk case, we examine here the influence of the number of scales retained for the construction of the frame and controller by b. For b = 5, 4, 3, the frame contains respectively 7, 8, 9 scales. The results are shown in Figure 2, where it can be seen that the frame composed of 7 scales produces the best results. Note that these performances might be improved, for example by making the SURE depend on the β parameter characterizing the threshold rule. The frame considered here has only one parameter, other more flexible constructions, based on a partition of the unit or other types of tight frames such as spectrum adapted and/or signal adapted tight frames from [START_REF] Shuman | Spectrum-adapted tight graph wavelet and vertex-frequency frames[END_REF] and [START_REF] Behjat | Signal-adapted tight frames on graphs[END_REF] could also lead to an improvement.

Pittsburgh Census Tract Graph

For the sake of completeness, our methodology is also compared to the trend filtering on the Pittsburgh Census Tract graph considered in [START_REF] Wang | Trend filtering on graphs[END_REF] which consists of 402 vertices and 2382 edges. The very same piecewise linear function ω with b = 2 is still used for this graph.

The number of scales is then J + 1 = 7. For this experiment, only the level dependent threshold procedure is considered.

We consider the signal and 10 realizations of the noisy signal generated in [START_REF] Wang | Trend filtering on graphs[END_REF] (corresponding to an average noise level of 4.84 ± 0.37dB). The resulting SNR for the oracles of SGWT and fused lasso (i.e. k = 0) are respectively 11.51 ± 0.52dB and 9.85 ± 0.54dB.

Additionally, we run a comparison with the trend filtering (i.e. for k = 0, 1, 2) for the signal f η,k with η = 0.01 and k = 5 with the different noise levels σ = 0.004, σ = 0.005 and σ = 0.01. A comparison with another wavelet estimator proposed in [START_REF] Sharpnack | Detecting activations over graphs using spanning tree wavelet bases[END_REF] is also provided, considering two thresholding rules (i.e. "soft" and "hard"). For these 5 competitors we only report the oracles results.

Even though the SURE no longer depends on the original signal, it does depend on σ 2 in both methodologies. Since in real applications, the noise level remains unknown in general, we introduce two naive estimators of σ. In fact, a straightforward computation shows that for any function g : R + → R + :

E[ f T g(L) f ] = f T g(L)f + E[ξ T g(L)ξ] = f T g(L)f + σ 2 Tr g(L),
so that a biased estimator of σ 2 is given by

σ2 1 = f T g(L) f Tr g(L) .
As soon as the original signal is reasonably smooth so that f T g(L)f is negligible compared to Tr g(L), then σ2 is an accurate enough estimation of σ 2 . As a first choice, we choose g(x) = x.

Thanks to Dirichlet's formula, it follows:

σ2 1 = f T L f Tr L = i,j∈V w ij | f (i) -f (j)| 2 2 Tr L .
This is nothing but the graph analogue of the Von Neumann estimator of von [START_REF] Von Neumann | Distribution of the ratio of the mean square successive difference to the variance[END_REF] explaining the terminology Graph Von Neumann estimator (GVN).

A second natural choice is given by g(x) = ψ J (x) corresponding to the filter at the finest scale. The resulting estimator is called High Pass Filter Von Neumann (HPFVN). The value of the estimator is easily computed from the coefficients as follows:

σ2 2 = n(J+1) i=nJ+1 (W f ) 2 i Tr ψ J (L) .
Table 2 summarizes the findings with the nomenclature of Table 1 (where LD stands for level-dependent (LD) coordinatewise thresholding) of the main document. Lines SURE β=1,2 σ1,2 stand for the SURE procedure in which the noise level is estimated by the GVN (σ 1 ) and the HPFVN (σ 2 ). HPFVN and GVN provide a very comparable sigma estimate in this setting.

Additionally, a visual comparison of our methodology with the fused lasso is illustrated in Figure 4 (corresponding to one realization in the context of the third column of the Table 2).

We can see that our approach provides a gain of about 2.5dB compared to the fused lasso. Again, in these experiments, the multiscale analysis shows better performances than the trend filtering. Besides, the estimation of σ is sufficiently accurate to have a fully data-driven procedure.

To conclude, let us stress that the SGWT is computationally more efficient than the GTF (k = 1, 2). For the latter, on a standard laptop (Intel Core i7@2.7GHz-16Go LP-DDR3@2133MHz), each 10 realizations consumed about 4h16m cpu time in mean (min:1h45m, max:6h19m) with the same grid search as for the Pittsburg in [START_REF] Wang | Trend filtering on graphs[END_REF]. The CPU consumption for k = 0 is more decent with a mean of 3m for each 10 realizations exploiting the idea and the c++ code of [START_REF] Chambolle | On total variation minimization and surface evolution using parametric maximum flows[END_REF]. Incidentally, this main drawback of GTF was noticed in [START_REF] Padilla | The dfs fused lasso: Linear-time denoising over general graphs[END_REF] forcing a preprocessing of the graph using Depth-first search algorithm.

The SGWT on its side consumed 3min for the diagonalization and 42s for each 10 realizations.

Real Dataset: New York City Taxis

Our methodology has been also tested on a real data fetched from NYC taxis 5 databases. We build a graph with 265 vertices consisting of the LocationID (Pick-Up and Drop-Off) and define Gaussian weights w ij = exp(-τ d 2 i,j ) where d i,j is the mean distance taken on all the trips between i and j or j and i. The signal f considered is defined upon the variable "total amount" on which an artificial noise is added. For an average input SNR in = 5.23 ± 0.38dB on N = 25 observations, we obtain for the level dependent SURE with β = 2 and σ known, an output SNR The quality of the SURE for uncorrelated noise is closely related to the intensity of the correlation tuned by the parameter α. As an example, if we choose α = 0.1, the SURE adapted to correlated noise still performs very well with an estimated SNR of 10.44dB compared to the oracle 10.46dB. The SURE for uncorrelated noise is nonetheless not that bad since it estimates the SNR at 9.78dB.

Consequently, the SURE estimate for uncorrelated noise is robust to small correlations which is particularly interesting in applications since it can be difficult to estimate the correlation structure.

Further Experiments with Block Thresholding

Finally, we report some experimental results in the context of block thresholding in the setting of Section 4.3. For each scale j = 0, 1, . . . , 7, the n wavelet coefficients are split into L blocks of uniform length (except for the last block that can be shorter). The best performance of block thresholding is achieved for blocks of size |L| = 47. The SNR in = 9.64 ± 0.33dB for 25 realizations. For the oracle global coordinate-wise threshold with β = 2, we obtain a SNR of 11.07 ± 0.33dB compared to the block procedure with a uniform threshold 11.82 ± 0.40dB.

The block procedure performs better than the coordinate-wise one for a uniform threshold but is actually worse compared to the level dependent coordinate-wise thresholding process.

The level-dependent method is expected to give better performance, but would require a more sophisticated optimization algorithm than grid search to be computationally acceptable.

Conclusion and Perspectives

In this paper, we have introduced a version of the SURE designed for SGWT, allowing automatic parameter selection in denoising tasks of signals on graphs. Closed-form expressions for coordinatewise and block SUREs have been provided for a wide range of threshold rules. Finally, the case of a correlated noise in the graph has also been considered. Many experiments on the Minnesota graph, the Facebook graph built, the Pittsburgh graph and the NYC taxis graph from real data are conducted.

For signals of different regularity on those graphs, it has been shown that the SURE provides an efficient estimate of the theoretical MSE. These experiments also show that multi-scale analysis is a serious competitor to existing methods. Indeed, the SGWT shows performances equivalent and sometimes even much better than the GTF, especially at high noise levels.

To be complete, the variance parameter σ 2 should be estimated. This has been (partially) addressed in the supplementary material by introducing the two estimators HPFVN and GVN.

  section presents the empirical performance of the proposed automatic threshold selection for different signals defined on different graphs: the Minnesota roads graph (seen as a reference in many recent studies, see Behjat et al. (2016) and references therein) with synthetic signals and the Facebook graph with signals from Wang et al. (2016), the Pittsburgh Census Tract graph, a graph built from a dataset on New York City taxis with a real signal as well as numerical experiments in the correlated and block cases. All the experiments are conducted with the R package gasper.

Figure 1 :

 1 Figure 1: Signals used on the Minnesota graph.

  p.14 and Figure 9 p.27)). With the exception of the dense Poisson case, where all methods provide Noise Level: -SnR in (dB)

Figure 2 :

 2 Figure 2: Mean SNR performance over N = 5 realizations on the Facebook graph.

4Figure 3 :

 3 Figure 3: Influence of the number of scales on the denoising performance.

Table 1 :

 1 Mean SNR performance over N = 10 realizations of the low to high noise levels settings with corresponding empirical standard deviation. Left panel: f 0.01,2 and right panel: f 0.001,4 .

	SNR in 16.07±0.13 10.05±0.13 4.03±0.13 16.64±0.13 10.62±0.13 4.60±0.13
	MSE β=1,G 19.04±0.24 14.22±0.26 9.46±0.26 24.67±0.33 19.77±0.36 14.79±0.45
	MSE β=2,G 20.07±0.24 15.60±0.30 10.69±0.30 26.88±0.29 22.18±0.37 17.08±0.67
	SURE β=1,G σ	18.96±0.27 14.16±0.29 9.46±0.26 24.62±0.37 19.64±0.48 14.70±0.52
	SURE β=2,G σ	20.04±0.37 15.49±0.43 10.64±0.32 26.73±0.32 21.91±0.52 16.88±0.59
	MSE β=1,LD 19.10±0.24 14.28±0.27 9.58±0.27 24.68±0.34 19.79±0.36 14.83±0.46
	MSE β=2,LD 20.08±0.24 15.61±0.29 10.72±0.30 26.90±0.26 22.20±0.36 17.13±0.67
	SURE β=1,LD σ	19.10±0.24 14.26±0.26 9.48±0.24 24.51±0.40 19.59±0.46 14.69±0.49
	SURE β=2,LD σ	20.01±0.31 15.51±0.36 10.61±0.39 26.52±0.32 21.79±0.34 16.73±0.60
	Wiener 17.01±0.13 11.87±0.15 7.43±0.16 17.91±0.12 12.86±0.13 8.40±0.16
	r inf 17.05±0.00 11.89±0.00 7.42±0.00 17.96±0.00 12.91±0.00 8.46±0.00
	MSE k=2 17.35±0.13 11.43± 0.16 5.68±0.17 19.37±0.14 13.65±0.18 8.01±0.24
	MSE k=1 18.05±0.14 11.98±0.18 6.24±0.18 20.43±0.15 14.78±0.18 9.30±0.27
	MSE k=0 19.57±0.17 13.43±0.23 7.62±0.22 23.38±0.25 17.88±0.27 12.80±0.40
	Generally speaking, we observe from Table

Table 2 :

 2 Mean SNR performance over N = 10 realizations with f 0.01,5 for the Pittsburgh graph.

	σ	0.004	0.005	0.01
	σ1	0.0065	0.0072	0.0113
	σ2	0.0068	0.0074	0.0114
	SNR in	9.53±0.29 7.60±0.29 1.58±0.29
	MSE β=2,LD 13.65±0.29 12.22±0.34 8.46±0.41
	MSE k=2	12.28±0.25 11.05±0.22 8.14±0.29
	MSE k=1	13.10±0.21 11.78±0.22 8.09±0.36
	MSE k=0	12.83±0.23 11.42±0.24 7.56±0.46
	MSE Soft	11.22±0.23 9.62±0.27 5.22±0.22
	MSE Hard	10.22±0.28 8.50±0.38 4.03±0.14
	SURE β=2,LD σ	13.33±0.37 12.00±0.30 7.95±0.25
	SURE β=2,LD σ1	12.35±0.53 11.38±0.62 8.13±0.35
	SURE β=2,LD σ2	12.12±0.57 11.27±0.54 8.07±0.33

https://github.com/fabnavarro/gasper

https://github.com/fabnavarro/SGWT-SURE

Available here: https://sites.cs.ucsb.edu/ yuxiangw/resources.html

https://s3.amazonaws.com/nyc-tlc/trip+data/yellow_tripdata_2018-01.csv

Correlated Noise

On the Minnesota roads graph and for the signal f η,k of Section 4 with η = 0.01 and k = 2, we add one realization of a correlated noise with covariance matrix given by (3) where α = 0.5 is added leading to a noisy signal with SNR in = 2.07dB. We run the level-dependent coordinatewise thresholding process. The SNR given by the oracle involving the unknown signal f η,k and the one given by the SURE estimator adapted to correlated noise are similar: 8.54dB versus 8.51dB. In this case, the SURE for uncorrelated noise shows very bad performances: we found 4.33 for the corresponding SNR.

The performance of these estimators highly depends on the underlying signal. Further investigations seem to be necessary.

Theoretically, the coordinatewise soft thresholding in the transformed domain satisfies an oracle inequality as shown in [START_REF] Göbel | Construction of tight frames on graphs and application to denoising[END_REF]. Similarly to the regular case, this oracle inequality states that the estimator mimics the oracle projection up to a log factor. The proof relies on the fact that the multivariate risk is expressed as a sum of univariate risks so that the Donoho's machinery applies. Using this fact, a maximal inequality for the SURE might be stated as well.

Regarding, the numerical complexity, the main limitation is the need of a complete reduction of the Laplacian. In the same vein as [START_REF] Hammond | Wavelets on graphs via spectral graph theory[END_REF], many of the involved steps might be numerically optimized using Chebyshev polynomials. Actually, the only problematic step in the method is the computation of the weights (WW * ) i,j appearing in the SURE. However, their expression in terms of covariance suggests that Monte-Carlo estimation could work. Besides, the space-time complexity might be reduced taking advantage of the low-rank property of WW * implying several linear constraints on the weights (precisely nJ). Finally, for the block procedure to be completely useful, an adapted optimization algorithm should be implemented. of frame should involve a notion of graph limit such as the one introduced for graphons (see [START_REF] Lovász | Large networks and graph limits[END_REF]) or the more probabilistic Benjamin-Schramm limit introduced in [START_REF] Benjamini | Recurrence of distributional limits of finite planar graphs[END_REF]. This formal study should also give rise to less nave estimators of the noise level and above all give recommendations according to the class of signals considered.