Plasmonic-organic hybrid electro/optic Mach-Zehnder modulators: from waveguide to multiphysics modal-FDTD modeling

Plasmonic organic hybrid electro/optic modulators are among the most innovative light modulators fully compatible with the silicon photonics platform. In this context, modeling is instrumental to both computer-aided optimization and interpretation of experimental data. Due to the large computational resources required, modeling is usually limited to waveguide simulations. The first aim of this work to investigate an improved, physics-based description of the voltage-dependent electro/optic effect, leading to a multiphysics-augmented model of the modulator cross-section. Targeting the accuracy of full-wave, 3D modeling with moderate computational resources, the paper presents a novel mixed modal-FDTD simulation strategy that allows to drastically reduce the number and complexity of 3D-FDTD simulations needed to accurately evaluate the modulator response. This framework is demonstrated on a device inspired by the literature. © 2020 Optical Society of America under the terms of the OSA Open Access Publishing Agreement


Introduction
Short range, high-speed communications probably are the most critical issue for the formulation of the forthcoming ICT strategies. Indeed, the present zettabyte age has been possible thanks to the centralized computation paradigm based on data centers, which has been enabled by high-speed, low-power optical interconnects [1]. This is why short-range communications are becoming as important as their consolidated long-haul counterparts, still necessary to deliver services at the end users [2]. Silicon photonics (SiPh) is one of the most viable platforms towards the hypothetically-coming yottabyte age, with its promise of the synergistical, low-cost and CMOS-compatible integration of electrical and optical systems [3].
Remarkable examples of SiPh-compatible devices are plasmonic-organic hybrid (POH) electro/optic (E/O) devices such as Mach-Zehnder [4] or ring [5] modulators for 1.3 m and 1.55 m communications, or disk resonators [6] for optical neural networks for deep learning. Focusing on Mach-Zehnder modulators (MZMs), the non-diffraction limited characteristics of the plasmonic waveguide enable nanoscale cross-sections and microscale total lengths, allowing chip-scale integration [7,8]. Such small cross-sections lead to very large radiofrequency (RF) electric fields with reduced driving voltages, enhancing the E/O effect and allowing for sub-THz bandwidths. These exceptional features are paid with the very strong propagation losses characterizing plasmonic modes, which are about 1 dB/ m. Nevertheless, the extremely compact achievable footprints enable an energy consumption of the order of fJ/bit, making these devices attractive for low-power communication systems [9][10][11].
In POH modulators, the E/O material is based on chromophore molecules dispersed in a host polymer medium, that are previously oriented by a static poling electric field [12][13][14]. Modulation of the material refractive index is enabled by applying a RF electric field to the poled material. This material fills the phase shifters slots, which are designed to support plasmonic modes [15]. It is therefore clear how a comprehensive model should predict the electro-optic modulation from RF electrical simulations, whose results are used to obtain a complex, position-dependent refractive index profile as the input of the optical model. A simplified model based on this principle could exploit 2D modal simulations to obtain the modulator static/dynamic response [16]. Yet, such an approach cannot account some effects, such as slot coupling and power losses related to surface plasmons pertinent to the top/bottom waveguide walls. Describing these details would require a 3D full-wave model of the entire geometry, which is very challenging because of the extremely severe memory and computational requirements.
The aim of this work is to present a critical appraisal of the available simulation approaches for POH MZ modulators. In Section 2 a reference device inspired to the literature is described (with the geometry details reported in Appendix A) and waveguide-only simulations are presented for both the cold regime (no applied RF input voltage) and electro/optic operation, opening a discussion of advantages and shortcomings of semi-analytic approaches compared to more sophisticated mode simulations, and emphasizing the importance of a multiphysics simulation framework. The only way to overcome the limitations of waveguide-based models is to move towards 3D full-wave modeling, leading however to huge computational costs and times. To mitigate the issue of CPU intensity, Section 3 presents an efficient multiphysics 3D approach based on a mixed modal-FDTD (MFDTD) simulation strategy, yielding accurate simulations with a drastic reduction of the computational burden, also made possible by the mode-matching technique formulated in Appendix B. Section 4 draws some final remarks on this work.

Multiphysics waveguide simulations
The device under analysis, sketched in Fig. 1, is a POH E/O MZM similar to the one presented in [19]. The structure is fabricated on a SiO 2 layer 3 m thick, grown on a Si substrate (not shown in the figure but included in the electrical simulations). The two arms of the MZM are the slot waveguides embedded between the central gold island and the two lateral gold rails. The optical input signal is assumed to be the fundamental mode of the input (left) Si waveguide. This mode reaches a splitter, consisting of a couple of facing tapers (left in Si, right in Au), which convert the dielectric waveguide mode into the plasmonic modes supported by the slots. The device is symmetrical with respect to the central section (the center of the gold island), so that, after propagating in the slots, the plasmonic modes are recombined and couple to the The blue curve has been simulated with an electromagnetic mode solver based on FEM [17], including all the geometrical details presented in Fig. 1(right) for ≥ 0. The red curve has been obtained approximating the geometry as a metal-insulator-metal waveguide and using semi-analytical expressions from [18,Ch. 10].
output Si waveguide. The device is immersed in the DLD-164 non-linear optic (NLO) material, with thickness ℎ NLO . Modulation is achieved through the electro/optic effect induced by the RF voltage, which is applied to the central island contact. The RF field changes the effective refractive index of the plasmonic modes in the phase shifters, leading to a voltage-dependent interference at the output combiner, which ranges from constructive (ON state) to destructive (OFF state).
The phase modulators are driven in push-pull operation by a single signal, using the coplanar ground-signal-ground transmission line sketched in Fig. 1(right) [8,[20][21][22]. This is obtained by aligning the poling field for the E/O polymer to the modulation RF field, the latter having opposite polarity in each modulator arm. By inspecting the geometrical details reported in Table  1, one can notice that the two slot widths are different. This comes from a precise choice, as it allows to tune the MZM to operate around the quadrature point (where linearity is maximum) at zero bias voltage [19].

Cold (zero voltage) device simulations
A preliminary investigation of these aspects can be performed on the basis of Fig. 2(top), which shows the effective index of a slot waveguide, at zero applied voltage, as a function of the slot width. Here, the results of two models are presented at = 1.55 m. The blue curve has been obtained by simulating the 2D cross-section of an isolated slot (i.e., simulating the ( , ) cross-section shown in Fig. 1(right), just for ≥ 0). The simulation of plasmonic slot waveguides has been widely addressed in the literature, with a broad spectrum of techniques including the effective index method [23], circuit approaches [24,25], finite difference schemes either in time [26] or in frequency [27] domains, finite elements [28], Fourier modal methods [29], and integral-equation schemes [30]. In this work, waveguide simulations have been performed img\src\neff_vs_VRF_and_SysLevel\05 -2020-05-28 -Tibaldi risposta ES\Main_0 vari Simulazioni non solide: NO multifisico ( Primo scopo figura: mostrare impatto d slot: le curve dashed sono ottenute sim mentre quelle dash-dotted tengono con modale, simulando entrambe le guide, di anticrossing, tipico nella coupled mo le guide sono vicine; qui, è più un mode diverse guide. Curve SOLIDE: il modello RF è valutato c (cit. Bertazzi). Notare che a 0 V le curve giustamente coincidenti: no E/O e quin with an in-house electromagnetic mode solver based on the finite element method (FEM) [17]. On the other hand, the red curve has been obtained with a much simpler and widely-available model, based on approximating the slot geometry as a metal-insulator-metal waveguide (therefore, -invariant) and using semi-analytical expressions from [18,Ch. 10]. From the top panel one could deduce that the simpler model, even though capturing the general trend vs. the slot width, is inadequate for design purposes. This is partially contradicted by the analysis in the bottom panel.
Here, recalling that the phase shift of a single MZM arm at 0 V applied voltage is Φ = 0 eff mod , we plot the phase difference between the two arms (where one of the arms is 100 nm wide), as a function of the width of the other arm: where 0 = 2 / , and the phase modulator lengths are assumed to be mod = 6 m. Obviously, ΔΦ = 0 • at slot = 100 nm, which corresponds to the symmetric MZM case. While the semi-analytic model fails to predict ΔΦ for wider slots, it is pretty accurate for narrower cases. Targeting at ΔΦ = 90 • , i.e., setting the half-power point at RF = 0 V, both models predict slot 85 nm (indicated with the blue and red open bullets), with a deviation smaller than 1 nm. Intuitively, this partial success of the semi-analytic approach could be ascribed to the slot aspect ratio (ℎ slot = 220 nm). Being the slots narrow, they are quite similar to metal-insulator-metal (MIM) waveguides, justifying the partial validity of the simplification that results useful to perform preliminary optical characterizations of passive plasmonic slots before their electro-optic implementation in a POH MZM [31].

Multiphysics-augmented waveguide simulations
The results discussed in the previous subsection pertain cold device operation; because of the absence of the modulating radiofrequency field, optical-only simulations are sufficient in this case. Modulation is enabled by imposing an RF electric field, which causes the optical dielectric where NLO is the element of the permittivity matrix relating the optical electric field along the poling direction with the displacement field along the same direction. To achieve (2), we neglected the -dependence of the electro-optic coefficient. This is acceptable in the splitter/recombiner, where the electric field profile is much weaker than in the slots. It is to be remarked that, even if the device under study is simplified (vertical slot walls, isotropic permittivity), in the transverse ( , ) plane this model can describe complex geometries, e.g., including slanted walls such as in [32, Fig. 10(c)], and sophisticated electromagnetic properties, e.g., position-dependent anisotropic permittivity. On the other hand, (2) ignores a possible -dependence of the poling field, which could arise for example from fluctuations of the slot width/height. The experimental characterizations (see the SEM image from [32, Fig. 10(b)]) suggest that it is reasonable to simulate an average slot width. With this hypothesis, the device geometry can be designed on the basis of parametric simulation campaigns, and only the final design verification could be performed by a full 3D electro-opto simulation, limiting the overall computational cost. Because the modulator length is much shorter than the RF wavelength, its frequency response can be reasonably approximated with that of a circuit, being the device and driver total equivalent resistance and its static capacitance. This allows to reduce the electrical analysis to a quasi-static problem in the 2D cross-section, as shown in Fig. 1(right), and to introduce this -independent E/O effect only in the phase modulators. Here, Δ mat can be evaluated as [33]: where 33 is the component of the E/O tensor (in contracted index notation) that affects mat due to an RF electric field applied along the poling direction (consider that in a quasi-TEM approximation the RF and DC field patterns coincide), and sign ( ) takes into account that it has opposite signs in the two slots (refer to the RF circuit sketched in Fig. 1(right)). It is understood how (3) leads to a multiphysics treatment, since the position-dependent ,RF , ,RF field components should be assessed through electrical simulations. Then, they should be interpolated on the optical problem mesh to evaluate Δ mat in its cross-section (the one simulated to produce Fig. 2), thus requiring a coupled, multiphysics approach.

V
To avoid this sort of complications, recalling that in MIM waveguides plasmonic modes are TM, with a dominant transverse component of the optical field, and that the RF field in the slots is mainly orthogonal to the slot walls, i.e., only the perpendicular ( ) component survives, the following approximation is sometimes adopted (see, e.g., [31]): yielding a constant, non-zero Δ mat only in the slot. Figure 3 shows the effective refractive index versus the RF voltage RF . There are three groups of curves, obtained with different degrees of approximations. The green dashed curves are obtained simulating the two slot waveguides separately, treating them as isolated just like in Fig. 2, and the approximated E/O refractive index variation model in (4) has been adopted. These curves clearly are straight lines, intersecting at RF 4 V, which corresponds to the MZM ON state (i.e., the two slot line optical fields are in phase, having the same optical path length).
a cross-section including both slots. Rather than being straight lines, these curves exhibit an almost-parabolic behavior in the proximity of the ON state, and become linear far from it, with almost the same slope of the dashed green lines. Because the two eff curves are not intersecting, the corresponding phase shift is not zero, suggesting that operation at RF 4 V should be a quasi-ON state, characterized by an excess loss; this point is further discussed in Section 3.3. Finally, the solid blue curves have been obtained including mode coupling and the full description of the E/O effect from (3). Here, the RF fields have been simulated with our in-house quasi-static (QS) FEM electric solver [34]. Compared to the dash-dotted red curve, the multiphysics simulation result exhibits a higher slope, which makes the ON voltage to be shifted at about 3.5 V. In this view, it is clear that (4) underestimates the E/O effect. The parabolic behaviour of the red and blue curves is caused by coupling effects between the two slot modes. Despite in the device under study the slots are quite far away, separated by the gold island, mode coupling is fostered by the surface plasmons of its top/bottom walls. Such behaviour, commonly referred to as anticrossing, is indicative of coupling between two modes [35]. This occurs in several EM structures, such as dielectric and photonic crystal waveguides [36,37] and high-contrast gratings [38][39][40], but also in semiconductor crystals where coupling between waveparticle modes is present [41].
The strong losses of the top surface plasmonic mode make the anticrossing strongly dependent on the slot separation island . This is investigated in Fig. 4. In particular, the left panel shows the eff ( RF ) characteristics curves for three different island widths. Here, the blue curve corresponds to the reference ( island = 400 nm) case reported with the same color in Fig. 3. Instead, the red curve ( island = 200 nm) exhibits a much broader eff splitting, Δ eff (definition in the figure), as a consequence of the increased mode coupling. On the other hand, the orange curves ( island = 600 nm) are almost intersecting, therefore tending to the isolated slot case of Fig. 3. A more quantitative information is provided by the right panel, showing Δ eff as a function of island . Waveguide simulations can be tentatively used to perform preliminary estimates of the MZM performance. To this aim, in Fig. 5, the electro/opto simulation results from Fig. 3 have been used to obtain a first estimate of the modulator . To this aim, the top panel shows the phase shift ΔΦ( RF ) between the two MZM modes at the recombiner section computed as where eff,1 , eff,2 come from the effective index curves reported in Fig. 3. The has been approximated as the difference between the OFF voltage, at which ΔΦ( OFF ) = 180 • , and the ON voltage, at which ΔΦ( ON ) is minimum (zero, in the isolated slots case). Remarkably, the estimate obtained with the multiphysics-augmented simulation (blue curve) are very close to the experimental findings discussed in [19,31].
To further clarify the importance of a multiphysics treatment, the bottom panel of Fig. 5 shows as a function of the modulator length mod , evaluated as in the top panel. From these results it is not possible to appreciate significant differences between the isolated (green dashed) and coupled (red dash-dotted) slot cases. Instead, it appears that the multiphysics simulation might lead to relevant variations of , which can reach almost 4 V for short modulators. These differences can be understood by analyzing the results reported in Fig. 6(left), showing the magnitude of the electric field simulated the the QS solver.
The slots can be identified as the regions where the electric field is stronger (tending to red). Moreover, it could be noticed that the field level is slightly higher in the left slot (i.e., the narrower one). The electric field is non-vanishing also out of the slots, which cannot be taken into account by (4). This can be better appreciated in Fig. 6(right), which shows field cuts performed in the slot center (blue, = 110 nm), at the slot top (red, = 220 nm), and the end of the NLO material (orange, = 300 nm). From the blue curves it can be seen that the estimate (4) is very accurate inside the slot (considering RF = 1 V, the left slot field is 11.1 V/cm, the right one is 10 V/cm). Looking at the red curve, one could notice that the optical field diverges at the slot corners [42,43]. However, the most relevant effect in this context is the non-vanishing field corresponding to the island, which is related to a residual field component associated to the island surface plasmons. This is even more evident at the end of the NLO material.

Efficient comprehensive 3D simulation
In the previous section, the voltage has been estimated only on the basis of the effective refractive indexes obtained from waveguide simulations, emphasizing how a multiphysics-augmented framework is instrumental to reproduce the experimental findings. Simple system-level models for the slot optical fields interference at the output combiners, also accounting for the plasmonic loss in each slot (see, e.g., [16,Sect. 6.4]), are customarily exploited to provide an estimate of other important device figures of merit, as the modulator insertion loss (IL) and extinction ratio (ER). However, such simple models neglect a number of effects related to both slot mode coupling and the detailed description of the splitter and recombiner sections, which feed the plasmonic phase modulators and extract the signal from it.
From this viewpoint, the maximum realism is provided by 3D full-wave simulations of the entire device, which in principle can be carried out by commercially-available electromagnetic simulators implementing the finite-difference time-domain method (FDTD), such as RSoft FullWave [44] and Lumerical FDTD Solutions [45] (all the 3D-FDTD simulations used in this work have been performed with the latter). Lumerical, starting from a defined input field source (in this case the Si waveguide mode in the input splitter), returns the position-resolved 3D profile of the vector electromagnetic field on the entire device. The MZM response can be defined by post-processing this information. Because the modulator is embedded within a complex optical system featuring grating couplers and other components which can filter out spurious contributions, it is reasonable to base this definition on the fundamental mode of the output waveguide. By projecting the total (3D) field on it, a mode transmission coefficient 21 can be defined, whose absolute value squared can be interpreted as a out / in .  Fig. 1, and the orange walls indicate the end of the simulation domain; below it, the S splitter scattering matrix is sketched (both figures present a dashed red line, indicating the splitter end; at its right a transmission line section with length port is included to enhance the simulation stability, to be de-embedded with (7)).

Introducing mixed modal-FDTD simulation and cold splitter characterization
Even though the compact footprint of POH MZMs makes such an all-in-one approach not impossible, it is rather prohibitive, considering that each RF voltage level requires a different 3D-FDTD simulation. By inspecting the top view of the device in Fig. 7, one can imagine an intermediate strategy between all-in-one 3D-FDTD and waveguide-only simulations: instead of modeling numerically the entire device, it is possible to focus the 3D-FDTD simulations only on the splitter/recombiner (which are actually equal, just mirrored), and the central part by means of waveguide simulations and transmission line theory. Since the phase modulators consist of two waveguides, this device resembles a bimodal Fabry-Pérot interferometer (BFPI) [46], whose concept is sketched in the bottom part of Fig. 7. Here, the cavity transmission lines describe the plasmonic modes supported by the MZM, whose voltage-dependent dispersion properties have been presented in Section 2 and now are be re-used to avoid to simulate by 3D-FDTD the phase modulators, minimizing the simulation number and cost, and are coupled mutually and to the outer ports (fundamental Si waveguide modes) by the scattering matrices describing the splitter/recombiner. Then, the device response can be computed as the cascade of the splitter, cavity and recombiner transmission matrices. The first step of this mixed modal-FDTD (MFDTD) strategy requires simulating, with the 3D-FDTD, only the section sketched in Fig. 7(right), at RF = 0 V. This starts from the input waveguide, and is terminated after a length port from the splitter end, as indicated by the dashed vertical line; this transmission line segment is introduced just to avoid problems that originate from terminating the device too close to metal corners. This simulation provides by itself interesting data, such as an estimate of the coupling losses due to photonic-plasmonic interference (PPI) (see [31, Fig. 11]). In the following sections it will be discussed how, combining this simulation to the waveguide-only analysis presented in Section 2, it is possible to perform a full-device simulation, similar to that presented in [31, Fig. 8(b)], but including also 3D effects such as mode coupling through the surface plasmons of the top island walls.
Considering only 2 internal modes to describe the modulator response, Lumerical allows to compute, through the S-parameter sweep functionality, a 3 × 3 matrix where, e.g., port 1 indicates the fundamental mode of the Si waveguide, and ports 2 and 3 the two plasmonic modes considered for the slot waveguides. This matrix can be re-arranged as Coherently with Fig. 7(left), the subscripts "o" and "i" are used to indicate the ports located outside and inside the modulator cavity, respectively. In (6), the superscripts "L" are used to remark that the scattering matrix blocks include the transmission line segment long port indicated in Fig. 7(right). In order to obtain the final 0 volt splitter matrix S 0 , one should de-embed such transmission lines, requires defining (at RF = 0 V) the phase shift matrix E port as: where eff,i , can be obtained from (possibly multiphysics) waveguide simulations analogous to those shown in Fig. 3. Finally, de-embedding is performed by applying:

Voltage-dependent mode coupling effects
The superscript "0" in (7) indicates that the matrices (7) are computed and valid only for RF = 0 V. Even under the hypothesis (stated at the beginning of Section 2) of introducing the E/O effect in the phase modulators only, the splitter/recombiner scattering matrix depends on voltage, since the ports and transmission line parameters are defined starting from modal basis of the phase modulators, whose elements are voltage-dependent. This is shown in Fig. 8, which reports the cuts of the real part of in the slot center ( = 110 nm) for the two modes. Invoking the molecular orbital taxonomy, the coupled slots supermodes 1 and 2 are antibonding-and bonding-like, respectively. Three RF values are considered: 0 V (green dash-dotted curves: cold regime), ON = 3.5 V (blue dashed curves), and 8 V (red solid curves). The mode topographies are clearly voltage-dependent. In fact, at RF = 0 V modes are mostly localized in the left and right slots, respectively. At RF = 8 V, modes are still localized, but with switched order. Instead, at ON modes tend to assume odd and even parities. This is a signature of mode coupling, just like the parabolic-like behaviour of eff ( RF ) in Fig. 3: coupling is strongest at ON and decreases at smaller or lower applied voltages. For this reason, the matrix S 0 from (7) differs from the S ii indicated in Fig. 7, as it does not take into account such mode coupling. Therefore, obtaining the voltage-dependent S from S 0 still requires a change of basis matrix W( RF ) from the modes computed for the cold device to those with a non-zero E/O effect: Computing the 2 × 2 change of basis matrix W( RF ) in principle requires a mode-matching technique, where the field continuity at the interface between cold and biased waveguide sections should be performed including the complete mode spectrum. This is the case occurring, e.g., in high-contrast gratings: even if their operation can be described just by 3 × 3 scattering matrices, computing their entries requires a mode-matching with a large number of modes to expand/project the fields at the bar-air discontinuities [47]. In the case of dielectric waveguides this is even more troublesome, since the electromagnetic problem is theoretically unbounded, and one should include also the continuum part of the mode spectrum. A possibility is mimicking free-space by closing the problem within a very large computational box, whose modal expansion is known analytically [48]. A handier, still very general approach is through 3D-FDTD simulations of a discontinuity between a plasmonic slot without/with E/O effect. In practice, this could be achieved by simulating two short transmission line segments, and then de-embedding the lines in a similar fashion to (7).
As it could be appreciated in the cuts reported in Fig. 8, the mode topographies are very localized, which is a signature of their plasmonic character. Moreover, voltage introduces just a mild dielectric discontinuity within the slots, without introducing any other significant change in the geometry. This explains why the mode topographies without/with voltage are so similar. In this view, in these devices, one can approximate the mode basis at a given voltage as a linear combination of the elements of the basis at a different voltage. This consideration allows the matrix W( RF ) to be derived from the coefficients of this linear combination, which can be obtained solving a least-squares problem, whose formulation is reported in Appendix B. It has been verified that the change-of-basis matrices obtained with the two methods agree well, enabling to use both approaches with similar results. The effect of the voltage-dependent mode coupling on the plasmonic modulator feed can be appreciated in Fig. 9. Here, the dashed curves are obtained plotting S oi ( RF ), which corresponds to the case of a very long modulator ( mod ≥ 40 m). The solid curves are obtained plotting (I − S ii S ii ) −1 S io for mod = 6 m, i.e., including cavity effects (which are instead negligible in long interferometers due to the high plasmonic losses). In both cases, the curves can be interpreted as the excitation coefficients of the modes (being, more specifically, the magnitudes of the progressive waves) at the onset of the modulator section, with/without cavity effects. It should be noticed how two groups of curves exhibit similar trends for every RF , and are almost equal in the ON region. This suggests that in opposite to OFF , which is strongly affected by the modulator length, the ON and the corresponding mode excitation coefficients are quite independent on it.

Evaluating the modulator response
Once S is determined, it can be used to find the matrix S as: where mod depends both on RF and on the modulator length mod : img\src\MZM_AllInOne\Main_MZM_response.m img\src\neff_vs_VRF_and_SysLevel\05 -2020-05-28 -Tibaldi risposta ESMain_2_2_Plot_Response Risultati: system leve impera. Nessun risca sono i risultati ottenu Dal momento che qu da modo fondamenta a guida dielettrica us ragionevolmente rap succede in un sistem dentro questo modul senso, si ha un IL di 4 come numeri sono pi Finally, the modulator reflection ( 11 ) and transmission ( 21 ) coefficients can be obtained by cascading the two matrices from (8) and (9): The validation of the MFDTD is performed versus the all-in-one 3D-FDTD results, as shown in Fig. 10. Each of the all-in-one simulations (one different 3D-FDTD simulation for each RF voltage) consist of 144 millions Yee nodes, requiring about 7 hours on a HP ProLiant DL560 Gen9 computer (featuring 512 GB RAM), parallelizing on all the Intel Xeon E5-4627 v3 (10-core) four CPUs. On the other hand, MFDTD requires a single 3D-FDTD simulation performed at RF = 0 V, involving just the splitter section (48 millions Yee nodes), which are combined to the multiphysics-augmented waveguide simulations to trace the full (voltage-dependent) modulator response.
The responses simulated with the two methods are reported with solid red and dashed black curves. The 3D-FDTD simulations (both all-in-one and splitter-only) have been performed using a uniform mesh (5 nm step in all directions) within the modulator section (central island, slots, part of the rails and gold tapers), and with the Lumerical auto non-uniform setting, with mesh accuracy parameter set to 5 (high accuracy). In the phase modulators, the E/O effect is evaluated from 2D quasi-static analyses (multiphysics approach) just like in Fig. 3, blue curve. The figure shows also the definitions of extinction ratio (ER) and insertion loss (IL), which are about 14.5 dB and 4.7 dB, agreeing qualitatively with the experimental findings from [19,31] (measured ER is 20 dB, measured IL is 5 dB).
The remarkable agreement between the two curves, even at −18 dB levels is validating not just the MFDTD algorithm, but demonstrates also the bimodal character of the MZM under study. It is to be remarked that the MFDTD algorithm could be extended, in a straightforward fashion, to  Fig. 11. Left: plot of the MZM response for different modulator lengths: the blue, red, and orange curves refer to mod = 3 m, 5 m, and 7 m, respectively. Right: plot of the extinction ratio (ER, thinner blue curve, referred to the left axis) and insertion loss (IL, thinner red curve, referred to the right axis) as a function of mod , obtained from 21 from the simulated modulator response (10). The thicker curves are used to emphasize the average trends. The top abscissas axis shows the corresponding .
devices whose operation involves a higher number of modes. In such cases, the rightmost matrix in (6) would be larger than 3 × 3, but it could be still possible grouping the parameters as L oo , S L oi , S L io , and S L ii . Even though these blocks would have different dimensions, eqs. (7)-(10) could be still applicable.
This method allows parametric investigations versus the modulator length at no additional computational cost. As an example, Fig. 11(left) shows the modulator responses obtained for mod = 3 m, 5 m, and 7 m (blue, red, and orange curves, respectively). Considering the longer modulators (smaller ), it appears that the curves are almost periodic. It is to be remarked that neither the orange nor the purple curves exhibit appreciable differences between the ILs evaluated at the two ON states, which instead could be expected from the waveguide simulations reported in Fig. 3. This could be understood from the analysis of Fig. 9: at the first ON voltage (about 3.5 V) the first, antibonding-like supermode, which in this case is quasi-odd (see Fig.8, top), is weakly excited. On the other hand, the bonding-like supermode is strongly excited. Being one supermode almost suppressed, no interference between the supermodes takes place at the recombiner, and the field propagates in both slots only according to the lower eff from Fig.3. This is why no excess loss can be observed. At the higher ON , the modes are strongly decoupled, i.e., localized in the single slots (as one can infer from Fig.8), and the output power results from their constructive interference at the recombiner, which in this situation can be estimated with (5).
A synthetic representation of the relation between the three fundamental figures of merit is reported in Fig. 11(right): , ER, IL. Focusing on the abscissas, Fig. 5(bottom) suggests that it can be, in first approximation, interpreted as 1/ mod , as reported on the upper horizontal axis. (Indeed, according to a simplified modulator model, the product mod is constant). The thinner curves exhibit an oscillatory behaviour, which is particularly evident in the ER curve, whose calculation involves the logarithm (dB) of small quantities. These oscillations can be ascribed to cavity effects (also visible in the results of Fig. 9); as a matter of fact, their amplitude decreases at increasing mod (and therefore at increasing losses in the phase modulators). Focusing on the average trends (thicker curves), the lower plasmonic losses in shorter modulators lead to reduced insertion loss. The interpretation of the trend of the ER (generally decreasing with and correspondingly increasing with mod ) is less obvious. Longer modulators are characterized by lower , therefore closer to ON , which is independent of mod . In this case, the excitation of the even mode at increases (approaching the mode 2 peak at 3.5 V in Fig. 9), thus deteriorating field extinction.
This analysis provides some guidelines towards the design of these devices. In MZMs, mode coupling is a detrimental effect, impacting in particular on the ER. As suggested by Fig. 4, which characterizes mode coupling on the basis of Δ eff at ON , better performance could be achieved for large island . As an example, the device presented in [49], exhibiting ERs greater than 20 dB, falls in this situation (slot separation of about 100 m). In order to achieve good ERs without increasing the transverse footprint, Fig. 11 suggests to sacrifice and design short interferometers. In this way the resulting modulator will have good IL, due to the moderate plasmonic losses, and reduced mode coupling at . As a different route, one could investigate coupler modulators where, as opposite to MZMs, mode coupling is the enabling physical mechanism [50]. Future works will deal with a comparison of the performance achievable by optimum MZM and coupler modulator designs.

Conclusions
The aim of the paper is to present a hierarchy of models for plasmonic-organic-hybrid electro/optic Mach-Zehnder modulators, applied to the simulation of a device inspired by the literature. The importance of a multiphysics approach has been investigated post-processing the results of waveguide mode simulations. Augmenting these models with electrical simulations of the RF field enables accurate predictions of , which are sufficient to characterize the E/O performance of the modulator.
Estimating the insertion loss and extinction ratio requires moving towards computationallyintensive full-wave simulations of the entire device. The all-in-one simulations produce results compatible with the experimental findings. Further refinements of the model would require detailed and accurate descriptions of the device, as well as of the material parameters (in particular, metal), which are not available at present and out of the scope of this work.
Aiming to reduce the number of 3D-FDTD simulations and their cost, a mixed modal-FDTD strategy has been formulated, invoking the analogy of the POH MZM with a bimodal Fabry-Pérot interferometer. The remarkable agreement between 3D-FDTD and MFDTD results demonstrates the bimodal character of the investigated device. In this work, the MFDTD has been applied to a specific device, but it could describe any Mach-Zehnder modulator. More in general, the method can in principle be extended to simulate modulators based on plasmonic rings. In fact, ring modulators can be divided into blocks (directional couplers and a ring phase shifter consisting in a single plasmonic slot), to be characterized either by a scattering matrix (the couplers) or through "azimuthal" (rather than standard) transmission line theory.

A. Geometry and material parameters
This appendix provides details about the geometry and material parameter used in this work. More specifically, Tables 1-2 contain the parameters of the geometry sketched in Figs.1 and  1(right).
The refractive indexes adopted in the simulations, mostly coming from typical literature values, are listed in Table 3. In particular, Au in the electrical simulation is treated as an impedance boundary condition, with conductivity = 4.1×10 7 S/m. Even though not visible in Fig. 1(right), the Si substrate has been included in the electrical simulations. The complex dielectric constant used for Au has been taken from [19], valid at = 1.55 . Also the DLD-164 polymer optical refractive index has been taken from the same reference. However, to the best of our knowledge no information is provided about its radiofrequency response (all the details we have found come from [32, Fig. 4]), so we assumed NLO = 1.83, with no dielectric losses, also in the quasi-static RF problem. Table 1. Geometrical parameters of the cross-section shown in Fig. 1 Value, unit 90 nm 100 nm 220 nm 300 nm 400 nm 520 nm Table 2. Parameters of the geometry shown in Fig. 1(left) not reported in Table 1.

B. Voltage-dependent change of basis: least-squares formulation
The approach described in this appendix could be seen as a mode-matching technique, where only two modes are used to represent the transverse field at the discontinuity. In other words, the modes of a waveguide subjected to E/O effect, | 1 , | 2 , are expressed as a linear combination of the zero-voltage modes | 1 , | 2 (the situation at which the splitter 3D-FDTD is simulated): Because this equation involves four coefficients , their determination requires formulating a 4 × 4 linear system, which is obtained projecting these two equations on two functions. In standard mode-matching techniques great attention is put on the projectors' definitions. In this case, being the system quite small, it is sufficient that the projecting modes are independent. For this reason, we choose to project the equations on 1 | and 2 |, leading to 1 | 1 = 11 1 | 1 + 12 1 | 2 2 | 1 = 11 2 | 1 + 12 2 | 2 1 | 2 = 21 1 | 1 + 22 1 | 2 where the bra-ket notation indicates the dot product where, as in the standard mode-matching, is the ( , ) cross-section of the waveguide discontinuity, the star superscript indicates complex conjugation, and A t , B t are the transverse fields, which must be continuous to satisfy the boundary conditions of Maxwell's equations. It can be noticed that the first and last groups of two equations are independent, leading to two uncoupled linear systems: The solutions of the systems (13), (14) are the elements of the matrix W in (8).