Nonlinear modeling of waveguide photodetectors

: A new method of simulating photodiode nonlinearities is proposed. This model includes the effects of non-uniform absorption in three dimensions, self-heating, and is compatible with circuit components defined in the frequency domain, such as transmission lines. The saturated output power and third order output intercept points of two different waveguide photodiodes are simulated, with excellent agreement between measurement and theory. The technique is then used to provide guidance for the future design of linear waveguide-based photodetectors.


Introduction
Microwave nonlinearities, specifically distortion and output power, have been modeled extensively in the time domain in the past. Typically, one or more sinusoidal tones excite either a semiconductor drift-diffusion model [1][2][3][4] of a photodiode or a nonlinear circuit equivalent [5]. The photodetector output is calculated for an appropriate number of cycles, a Discrete Fourier Transform is applied, and the output power and output intercept point are calculated. These models have proven successful in identifying the primary causes of nonlinearity in high-power detectors, which vary by cross-section and with operation frequency. Time-domain modeling has some drawbacks, however. The number of time samples required increases as disparate time constants are introduced, as would be the case in a system where thermal effects or low-frequency noise spectra were important. Relative to frequency-domain modeling, the computational cost of adding even a linear component to the circuit is high. This can make including impedance matching networks and non-uniform illumination effects prohibitively time-consuming.
In this paper, we demonstrate the use of the Harmonic Balance Method (HBM) to examine microwave nonlinearities in photodiodes. At the component level, the primary advantage of the HBM is its ability to quickly expand a well-understood cross-section to an additional spatial dimension. This is especially useful in modeling waveguide photodiodes. Additional advantages include reduced computation time for sensitivity analysis and the ability to accurately model circuits with multiple disparate time constants without additional computational complexity. For devices where self-heating plays an important role, this is a significant advantage. At the circuit level, HBM is compatible with circuit components that are best described in the frequency domain, such as transmission lines, whereas time-domain solutions are not. This enables the designer to predict the behavior of one or many high-power photodiodes in a realistic microwave circuit. This ability would be valuable in many systems applications, such as photodiode-driven phased antenna arrays, photonic oscillators, and photonic analog-digital conversion [6]. In each case, the system is benefitted by highly linear photodiodes, and hence, integrating a circuit-based photodiode microwave non-linearity model into the overall system model enables more efficient and accurate system designs. HBM has been proposed as a viable simulation technique for photodiode nonlinearities in the past [7], but has not been used to optimize linear behavior.
To demonstrate the effectiveness of the harmonic balance method, we simulate the saturated output power and third order output intercept point (OIP3) of two different waveguide devices with previously well-quantified cross-section behavior. First, the output power from [8,9] is simulated. Excellent agreement between simulation and experimental results are achieved. The effects of different design parameters and operating conditions on output power are studied, and cross-section design and absorption profile are identified as the most important factors in determining performance. Second, the OIP3 from [10] is simulated. Good agreement between simulation and experiment is again obtained, and it is shown that although waveguide device performance can theoretically approach uniformly illuminated device performance, it may require careful optical designs to do this while maintaining good bandwidth and efficiency.

The linear model
We use Agilent ADS to implement nonlinear circuit models of the photodiodes studied. Figure 1 shows the generic detector circuit model used for all device cross-sections. To deal with non-uniform absorption in the propagation direction, the detector is divided into a series of unit cells, and a traveling-wave model [11] is used to model the optical and electrical connections between the cells. This is illustrated in Fig. 2. The use of an equivalent circuit model implies very little about the underlying semiconductor physics, and thus it is equally valid regardless of material system. The primary requirement for an equivalent circuit model to be valid is that the quasistatic assumption holds, i.e. that the value of the capacitance, responsivity, etc. at any given time is a function only of the terminal voltage and current at that same time [12]. The cross-section circuit model was originally proposed to incorporate transit-time effects into bandwidth simulations [13] and has been shown to accurately model linear behavior [5,13]. Here, we convert it into a nonlinear model by transforming linear elements into nonlinear ones. This limits the nonlinearities considered to those caused by changes in transit-time limited bandwidth, output impedance, and responsivity. These have previously been identified [2,5] as the dominant causes of nonlinearity for a wide variety of photodiodes, and are adequate to describe the detectors considered here. The leftmost section of the cross-section equivalent circuit deals with optical propagation and absorption. The center section deals with transit-time effects, and the far-right with diode impedance. The absorption transconductance g m1 , along with the associated voltage V opt , provides coupling between the optical input and the photodiode. V opt thus represents the input optical power. The second transconductanceg m2 , along with the associated voltage V τ , provides coupling between the transit time and diode impedance related circuit sections. The product of C τ and R τ is the transit-time pole. The diode capacitance, parallel resistance, and series resistance, are denoted C pd , R p , and R s , respectively. In the simulations undertaken here, the value of V opt was chosen such that the input power was equal to the microwave power available to a lossless photodiode driving a 50 Ω load. Preserving a meaningful relationship between V opt and the output power is not necessary and places restrictions on R in , g m1 , R τ , and g m2 . For example, for large R in and at DC, the output current (I) of the photodiode is Z 0· g m1· g m2· R τ ,where Z 0 is the characteristic impedance of the microwave source (usually 50 Ω). The product g m1· g m2 is therefore η/ (Z 0· R τ ), where η is the quantum efficiency of the detector at DC. Alternatively, the values of g m1 , R τ , and g m2 can be chosen arbitrarily, and then R in adjusted to obtain the correct quantum efficiency. This was the approach pursued in this work. It is worth noting that a detector circuit that is mathematically equivalent to this one but with fewer components can be drawn. However, using a simpler circuit would come at the expense of more complex expressions for nonlinear quantities. Figure 2 shows several cross-sections connected in parallel to form a single circuit model for a waveguide device. Each unit cell consists of the circuit shown in Fig. 1 (represented as Y PD in the second unit cell) and traveling-wave impedance Z PD .This formulation neglects carrier transport that may occur in the direction of optical propagation. This could be represented by additional connections between unit cells, but accurate results were obtained without doing this. The traveling-wave impedance can be found from analytical [14] or finiteelement solvers. For many devices, including the ones considered here, traveling wave effects are not very important, but replacing it with a short can result in numerical instability. Poor agreement with experimental results was obtained when using incorrect values for linear quantities, especially the diode capacitance and the series resistance. The optical propagation portion of the circuit (S in Fig. 2) can be represented by equivalent scattering parameters. It can be considered a direct representation of the photocurrent available in the optical signal. As such, it can be represented as a series of equivalent S-matrices with zero reflection (S 11 = S 22 = 0) and transmission with both loss and phase shift. The off-diagonal elements of the S-matrix are then exp((-jω/v g -Γα)z) where ω is the microwave frequency, v g is the group velocity of the mode, and Γ and α are the confinement factor and (power)absorption coefficient, respectively. Resistors to ground placed at the intersections between adjacent optical blocks can be used to represent excess modal loss. Any current that flows from one optical block but not into the next optical block reduces the efficiency, and is thus equivalent to excess loss. Lowering the value of R in has the same effect. Multiple optical paths operating in parallel can be used to model multimode effects. For short detectors, where there is still a significant amount of power in the optical signal at the end of the device, a terminating resistor (R term ) is necessary to avoid (nonphysical) microwave reflections from impacting the model.

The nonlinear model
To incorporate nonlinear effects into the linear model, coefficients such as C pd become functions of voltages and currents. The nonlinear functions that describe these coefficients can be measured for devices with uniform illumination (surface-normal devices), simulated in TCAD software, or in some cases described using analytical formulas. Measurement is the most accurate of these three, but is not always feasible and may give less insight into the internal workings of the device. Analytical formulas are generally easier to work with than TCAD models, but can be less accurate, especially for highly band-engineered structures. Here, we use an analytical formula to predict output power and experimental results from surface-normal devices to predict the OIP3 of waveguide devices with the same cross-section. Using the harmonic balance method imposes some restrictions on the kinds of functions that can be used to describe nonlinear elements. To avoid numerical instabilities, all functions must be both continuous and continuously differentiable. Furthermore, they should give physically possible results beyond the range of normal device operation. Thus the diode capacitance is probably better described by a polynomial fit than a function of the form 1/x, as a simple parallel-plate capacitance model would suggest.
In principle, a single function can be used to describe each nonlinear quantity regardless of the desired output of the simulation. In practice, this can lead to needlessly complex expressions and long computation times. Nonlinear capacitance, for example, plays an important role in determining OIP3 [5], but is not a dominant factor in saturation photocurrent at low frequencies. Thus it can be neglected entirely in some simulations. Similarly, the best functions to describe current-and voltage-dependent responsivity are not the same for both OIP3 and saturation photocurrent simulations. For a large-signal model used to predict the maximum output power of the photodiode, the output current can be assumed to be linearly proportional to the input power until the saturation photocurrent is reached. Drawing from expressions developed by the transistor simulation community, g m2 is given by functions of the form: Where J linear is the expected current density given the responsivity and k is a parameter used to change the shape of the curve. In terms of current, For values of k less than or equal to 1, the output current can never be equal to the maximum current (unless I max = 0) because tanh(1) < 1.In this work, k was chosen so that the output current would be equal to its maximum value when I linear = 0.99*I max , i.e when k = a tanh(0.99)/0.99. This improves accuracy near the onset of compression, but reduces it for heavily compressed photodiodes, as it allows the output current to exceed its maximum value when I linear >I max . Equation (1) is meant to model signal clipping, and is unlikely to accurately describe device behavior in the linear regime. Thus for OIP simulations, a polynomial fit will usually provide a better fit to the data.

Saturation power
To illustrate HBM's ability to deal with non-uniform illumination and thermal effects, we simulated the large-signal compression current of the PIN photodiode studied in [9]. It is a 7.4x500 μm waveguide detector with a 700 nm thick Germanium intrinsic region and a 4.4 GHz bandwidth limited by the RC time constant. The diode admittance can be inferred from [8] and simulation parameters are shown in Table 1. The optical power is assumed to decay exponentially along the length of the device with one of two characteristic lengths. Given the difference in responsivity between a 50 μm long device and a 500 μm one, the absorption length for one polarization state (pol.1) can be inferred to be around 38 μm (260 cm −1 ). The responsivity of the shorter device is polarization-dependent, and the absorption length for the lower-responsivity polarization state (pol. 2) is around 57 μm (175 cm −1 ). The basic nonlinear model consists of a series of5μm long sections. Since the operating frequency (1 GHz) is much lower than the −3dB frequency at the photocurrent levels of interest, the dependences of C τ and C pd on current and bias voltage are ignored and only the saturation current of the diode is considered. The maximum output current density of each section is ( ) where v eff = v n v p /(v n + v p ) and v n and v p are the saturated electron and hole velocities (both 6·10 6 cm/s in Ge [15]), ε i is the dielectric constant of the intrinsic region, W i is the intrinsic region thickness, V bias is the bias voltage, V pd is the voltage drop across the device due to the photocurrent and load impedance, and V th is the threshold voltage. This equation can easily be derived in the same way as the equation for Kirk current density in a bipolar transistor, for example as in [16], assuming the photogeneration rate is uniform and the background doping is negligible [17], and has appeared in high-power photodiode literature in the past [18]. As the current density in the intrinsic region increases, the charged carriers (holes and electrons) screen the electric field. Whereas under low-injection conditions the field is roughly constant, under high injection the field distribution becomes parabolic with the minimum occurring in the intrinsic region. The maximum current density is reached when the minimum of the electric field drops below the value necessary for holes and electrons to be able to maintain their saturation velocities. The threshold voltage is defined by the value of the minimum field. In germanium, the critical field is around 2 kV/cm [15]. This would correspond to a threshold voltage of 0.14 V in this device given the intrinsic region thickness of 700 nm. However, capacitance-voltage curves indicate that the intrinsic region is not fully depleted with less than 0.4 V applied bias, so this value is used for V th . Figure 3 shows the measured and simulated RF output power at 1 GHz as a function of photocurrent for pol. 2 (Γα = 175 cm −1 ). Fig. 4 shows the −1 dB compression currents and output powers for both polarization states. Good agreement (within 3 mA and 1.5 dBm for both polarizations) is achieved for the full range of operating voltages, with no fitting parameters used. As can be seen by comparing Fig. 3(a) to Fig. 3(b), the value of k in Eq. (1) was chosen to maximize accuracy near the onset of compression.  Heating reduces the compression current and output power in two ways. First, the absorption profile changes as a function of temperature. Since the direct band edge of Ge is close to 1550 nm, the absorption coefficient increases with temperature as the bandgap narrows. Second, as the temperature increases, the saturated carrier velocities decrease. This in turn decreases the maximum output current density. To take these effects into account, the absorption transconductance (g m1 ) and second transconductance (g m2 ) become functions of terminal voltage and output current. The temperature of each section is calculated separately, so that the carrier velocities and absorption coefficient are allowed to vary along the device.
The thermal impedance (Z t ) of each 5 μm long section is calculated to be 3000 K/W using a finite-element model that agrees well with measured data [19]. Changes in absorption due to self heating are modeled assuming a linear dependence of the absorption coefficient on temperature. The confinement factor can be a function of temperature in some photodiodes, but the variation was simulated to be very small for this device. The initial transconductance is scaled: Since α ~sqrt (E g ), To account for the change in propagation loss, the input impedance (R in ) is also scaled: Taking values for temperature dependent bandgap narrowing from [15], dEg/dT = −0.58 meV/K. The maximum current density decreases as the hole and electron velocities decrease. Taking values from [15], the decrease in v eff is approximately linear from 300 to 600 K and is 29 m/s/K. It is worth noting that poor agreement between measurement and simulation was obtained when thermal effects were ignored. We can then use the model to identify the most important design variables in determining the output power. To do this, we varied a single design variable (length, width, thermal impedance, absorption profile, intrinsic region thickness) while keeping all others the same as in the previous simulations. The effect of scaling on linear and nonlinear circuit elements from Fig. 1 was taken into account, but the thermal impedance was assumed to be independent of device dimensions. This was done because thermal effects can be mitigated substantially by better heat-sinking without significant effect on other device parameters [20].The maximum output current and capacitance of a given section were scaled linearly with area according to Eq. (3) and the expression for a parallel-plate capacitor. The contributions of the n-Ge contact, the p-region under the diode, and the p-Si contact to the device series resistance were not known, and these scale differently, so values were chosen such that the n-contact and p-spreading regions both made significant contributions and the total resistance added up to the values given in [8]. Figure 5 shows the effect of device length andwidth on output power. At the maximum output power, the front part of the diode will be compressed and a significant amount of useful photocurrent will be generated in the back part. Thus very short diodes have much lower maximum output power than longer ones. Increasing the length beyond 200 μm has minimal impact on the output power because most of the light is absorbed in the first 100 μm. For very long devices, the maximum output power decreases because the RC limit approaches the operating frequency. The width has a similar effect on output power. Very narrow devices are much worse than wider ones, but the width cannot be increased too much before the increasing capacitance begins to have an impact on performance. The data points shown in Fig. 5 are meant to provide qualitative confirmation of the simulation result; the experimental conditions differed slightly from the simulation. It is well-known that there are thermal limits to maximum output power. These come from two sources:(1) the maximum output power decreases as a function of temperature; and (2), there is a maximum temperature a device can reach before catastrophic failure [21]. Improved heat sinking has a significant impact on devices that are limited by their membership in the latter category, but minimal impact on devices in the former. Figure 6(a) shows the maximum output power of the device as the thermal impedance is decreased to 0 and increased to 10x its original value, for a fixed bias voltage of 6 V. Only the value of the thermal impedance, and not any device dimensions, is changed. Though a tenfold increase in thermal impedance would clearly be bad, it would not have as much of an impact as halving/doubling the width or decreasing the length below the characteristic absorption length. The results in [20,22] are consistent with this. At 5 V bias, devices with the same cross-section had 1dB compression powers around 23 dBm, whether or not they were flip-chip bonded to a heat sink. The maximum output power obtained in [20] was significantly higher than the power obtained in Assumptions made regarding the bias voltage at which the detector will be operated have a significant impact on cross-section design. Though Eq. (3) implies that the saturation current should increase monotonically with shrinking intrinsic regions, this is not the case for output power due to the voltage swing across the load. For a fixed bias voltage, there is some optimum design point where the load swing and field screening effects are balanced, as shown in the blue curve in Fig. 6(b). If we take into account the fact that devices with thicker intrinsic regions can be biased at higher voltages without going into avalanche mode, then the output power will increase with increasing i-region thickness until the transit-time pole approaches the operating frequency, resulting in a much broader optimum design point as shown in the green curve in Fig. 6(b). It has also been established that dilute absorption profiles are preferable to abrupt ones for power handling. This effect is very clear from Fig. 7, as the output power drops from 25 dBm to 7 dBm as the Γα product is increased from 1 cm −1 to 1500 cm −1 (this roughly corresponds to a confinement factor increase from 2% to 50%, though the value of α is a function of wavelength and growth conditions [23]. Figure 7(a) may be misleading as a design curve, however. The photodiode studied here was very long (500 μm), so all the devices considered have quantum efficiencies above 90%. In a practical system, there would be a trade-off between absorption profile, bandwidth, and quantum efficiency. This is shown in Fig. 7(b), where the design frequency is varied from 10 GHz to 40 GHz and the software is allowed to choose the width, length, and absorption profile in order to maximize the output power. For unlimited input power, the most dilute absorption profile allowed is the best one, regardless of design frequency. When there is a power budget in play, quantum efficiency and bandwidth concerns become dominant, and sharper absorption profiles are preferable.

Linearity
The harmonic balance method can be used to predict the third order output intercept point (OIP3). The detector from [10] [5]. Both surface-normal and waveguide devices with the same active areas were fabricated on the same chip and characterized together. The measurement result from the surface-normal device was used to construct a nonlinear model for the device crosssection. That model was then used in conjunction with an optical simulation in order to predict the OIP3 of the waveguide device. The cross-section model was very similar to the one presented in [5]: the responsivity and the admittance were the only nonlinear elements considered. A polynomial fit was used to express the responsivity and admittance in terms of terminal current and voltage, and coefficient values were chosen by fitting to the surfacenormal data. Table 2 shows the simulation parameters used for the model, where V τ is in volts and I is in milliamperes and refers to the current generated by the second transconductance element. Absorption-related parameters are not listed because no attempt was made to maintain a physically meaningful relationship between V opt and the photodiode output. A distributed model using the same coefficients (appropriately scaled) was used to predict the OIP3 of the waveguide device. The photodetector was divided into ten sections of equal length (10 μm). The absorption profile was assumed to have a single decay constant (Γα) of 335 cm −1 , which is consistent with BPM simulations. Figure 8(a) compares the simulation result and measurement at 10 mA of photocurrent, and good agreement is achieved over the whole frequency range. Figure 8(b) shows the simulated OIP3 at 7 GHz and quantum efficiency of the device as a function of characteristic absorption length. As the absorption profile becomes more uniform, the OIP3 increases and approaches the simulated value for the surface-normal detector at this frequency. However, the quantum efficiency at this point is very low, implying that more creative solutions, such as photodetector arrays [25] or a variable confinement factor [26] may be necessary to achieve surface-normal level performance from a waveguide-based device.   Figure 9 shows a collection of simulated absorption profiles and associated photocurrents. The three photodiodes considered have the same quantum efficiency (96%) and physical dimensions. The blue curve corresponds to a detector simulated above and has a predicted OIP3 of 21 dBm. The yellow curve is a linearly increasing Γα product and has an OIP3 of 25 dBm. The green curve increases approximately exponentially and achieves nearly uniform current density through the device. The OIP3 is 27.5 dBm, which is a little less than the 28 dBm the model predicts for the uniformly illuminated lumped-element device. The discrepancy is due in part to the simulated RF loss from the traveling wave electrode (0.25 dB from the front of the device to the probe pads), and in part to the illumination not being truly uniform in the simulation. Fig. 9. Simulated (a) absorption coefficients and (b) normalized photocurrent in the direction of propagation. The blue curves correspond to uniform Γα product, the yellow curves to linearly increasing Γα product, and the green curves to nearly uniform photocurrent density.

Conclusions
This simulation technique accurately predicts the linear and nonlinear behavior of photodiodes, enabling direct quantitative comparison of the effects of different design decisions on the figure of merit of interest. At the component level, this allows for simultaneous optimization for frequency response, output power, and OIP3.Using the harmonic balance method in a SPICE environment creates new potential for careful microwave design. Because HBM is compatible with transmission lines, matching networks for load impedance optimization can be designed and realistic interconnects between multiple devices can be simulated.