Accurate and Comprehensive Spectrum Characterization for Cavity-Enhanced Electro-Optic Comb Generators

Cavity-enhanced electro-optic comb generators (CEEOCGs) can provide optical frequency combs with excellent stability and configurability. The existing methods for CEEOCGs spectrum characterization, however, are based on approximations and have suffered from either iterative calculations or limited applicable conditions. In this paper, we show a spectrum characterization method by accumulating the optical electrical field with respect to the count of the round-trip propagation inside of CEEOCGs. The identity transformation and complete analysis of the intracavity phase delay were conducted to eliminate approximations and be applicable to arbitrary conditions, respectively. The calculation efficiency was improved by the noniterative matrix operations. Setting the maximum propagation count as 1000, the spectrum of the center ±300 comb modes can be characterized with merely the truncation error of floating-point numbers within 1.2 s. More importantly, the effects of all CEEOCG parameters were comprehensively characterized for the first time. Accordingly, not only the exact working condition of CEEOCG can be identified for further optimization, but also the power of each comb mode can be predicted accurately and efficiently for applications in optical communications and waveform synthesis.


Introduction
Optical frequency comb (OFC) is composed of a series of equally spaced and phase coherent frequency components [1,2]. Its unique property in the frequency and time domains brings revolutionary development in the fields of precision spectroscopy [3][4][5], optical communication [6,7], waveform synthesis [8][9][10], and precision metrology [11][12][13] etc. Compared to the OFC generation schemes based on mode-locked lasers [14][15][16] and micro-resonator lasers [17,18], the electro-optic modulators (EOMs)-based OFC generators have some unique advantages [19,20]. OFCs with a high repetition rate up to tens of GHz can be conveniently obtained with a robust and compact setup. Moreover, the central frequency and the repetition rate of the generated OFC can be configured freely and independently [21]. Such advantages make it a perfect multiwavelength laser source for applications in the fields of optical communication and optical arbitrary waveform generation [22].
Limited by the weak EOM interaction strength, the OFCs directly generated by a single EOM suffer from the narrow span of the comb spectrum. The cascade of multiple phase and amplitude EOMs can broaden the OFC spectrum up to tens of comb modes and a few nm of spanning, but no more [23]. Highly nonlinear fiber can be applied to broaden the comb spectrum as well [24]. However, the system structure and size must be increased by adding optical amplification and pulse shaping units to excite the highly nonlinear effect.
The shape and phase of the broadened comb spectrum is hard to predict or characterize as well [22].
An alternative wide spectrum EOM comb generation method uses an optical resonator to enhance the modulation process. This method was first proposed by T. Kobayashi et al. in 1972 [25]. Since 1993, a series of great studies on this cavity-enhanced electro-optic comb generator (CEEOCG) have been conducted by M. Kourogi, L. R. Brothers, A. S. Bell, J. Ye, and U. Sterr et al. [26][27][28][29][30][31][32][33][34]. As a new type of comb source, the comb spectrum characterization is the priority. In [26], M. Kourogi proposed the first CEEOCG spectrum characterization method by iteratively accumulating the inter-mode power coupling. The influence of material dispersion was investigated and compensated to achieve a broader spectrum range as well [31,32]. As the number of the comb modes and iterative processes are both theoretically infinite, however, even Kourogi himself admits that this calculation is too complicated [26]. To simplify the process, only the ±15 adjacent comb modes were iteratively calculated with matrix operations [30]. However, the iterative accumulation principle of this spectrum characterization method significantly magnified the errors from the limited comb modes and truncations during signal processing. The comb of a wider spectrum from CEEOCG still cannot be characterized accurately [35,36]. To avoid the complex iterative calculations, Kourogi proposed an exponential approximation method to characterize the spectrum of CEEOCGs [26]. However, the abundant mathematical approximations severely limit the accuracy of the spectrum characterization. Therefore, the existing CEEOCG spectrum characterization methods all suffer from insufficient accuracy due to either the iterative error accumulation or the abundant mathematical approximations.
At the same time, our previous research has proved that the accurate identification of the CEEOCG working condition is required for the fine adjustments of CEEOCG to achieve its optimized performance [37]. However, the CEEOCG comb spectrum has not been comprehensively characterized for the different working conditions of CEEOCG. Consequently, the lack of an accurate and comprehensive spectrum characterization method prohibits the wider application of CEEOCGs.
In this paper, we propose an accurate and comprehensive spectrum characterization method for CEEOCGs by accumulating the optical electrical field with respect to the count of the round-trip propagation inside of the CEEOCG cavity. Different from the existing methods, the proposed method is free from the iterative calculation of the power coupling among the generated comb modes. More importantly, there is no mathematical approximation introduced in the derivation. The influence of the limited count of the round-trip propagation was analyzed in detail and the accuracy of different methods was compared. Based on the proposed method, the influence of all the parameters was investigated independently and jointly. To our knowledge, it is the first comprehensive analysis of CEEOCG based on a highly accurate spectrum characterization method.

The Existing Methods for the Spectrum Characterization of CEEOCGS
A CEEOCG consists of an EOM inside a spatial linear cavity or an integrated ring cavity. As these two structures share the same essential principle, we mainly discuss the former type in this paper. A schematic of a typical spatial linear cavity CEEOCG is shown in Figure 1. When a single-wavelength seed laser is incident into the CEEOCG, it oscillates inside the cavity and passes through the EOM multiple times. Driven by a radio-frequency (RF) signal, the EOM introduces phase modulation sidebands to the seed laser. This effect of sideband generation is greatly enhanced by the repetitive beam propagations through the EOM and ensures OFC generation. As the comb modes come from the sidebands of the EOM phase modulation, the repetition rate of the generated comb is determined by the frequency of the EOM modulation signal. To characterize the comb spectrum, the power coupling effect among the comb modes was analyzed for the first time by M. Kourogi in [26]. As the total spectrum of the generated comb is composed of all the modulation sidebands, the electrical field of the transmission output from CEEOCG can be expressed as: where Ek is the electrical field amplitude of the k-th order sideband, ν0 is the frequency of the incident seed laser and ωm is the angular frequency of the EOM phase modulation. When the beam propagates inside of the linear cavity for a complete round-trip, the power coupling among the sidebands can be expressed as: where t1 is the electrical field amplitude transmission coefficient of the input mirror, Ekin is the electrical field of the k-th order sideband contained in the input laser spectrum, r1 and r2 stand for the electrical field amplitude reflection coefficient of the front input and rear output cavity mirrors, respectively and ϕk is the round-trip phase delay of the k-th order sideband. For a sideband of the q-th order, the power coupling efficiency is described by the first kind Bessel function Jk−q(2β) of the (k−q)-th order, where β is the phase modulation index of the EOM. As there are two processes of phase modulation in a roundtrip transmission inside of the linear cavity, the modulation index in Equation (2) is doubled as 2β. For an integrated ring cavity with a single section of phase modulation, the modulation index should be set as β [35,36].
It should be noticed that all the generated sidebands contribute to the k-th order sideband after a round-trip propagation. Therefore, the accumulation of all the coupled power from each sideband is required in Equation (2). At the same time, the power of the other sidebands is simultaneously changed together with the k-th order. To analyze the generated comb spectrum precisely, therefore, an iteration calculation of the inter-mode power coupling has to be applied. In fact, even the proposer M. Kourogi himself admitted that such an iteration calculation is too complicated [26].
In [30], an approximation calculation was made with less than ±15 adjacent sidebands, instead of all the sidebands ideally. The solution of a 1024 × 1024 sparse matrix equation was applied to describe the CEEOCG output spectrum. However, there are two To characterize the comb spectrum, the power coupling effect among the comb modes was analyzed for the first time by M. Kourogi in [26]. As the total spectrum of the generated comb is composed of all the modulation sidebands, the electrical field of the transmission output from CEEOCG can be expressed as: where E k is the electrical field amplitude of the k-th order sideband, ν 0 is the frequency of the incident seed laser and ω m is the angular frequency of the EOM phase modulation. When the beam propagates inside of the linear cavity for a complete round-trip, the power coupling among the sidebands can be expressed as: where t 1 is the electrical field amplitude transmission coefficient of the input mirror, E kin is the electrical field of the k-th order sideband contained in the input laser spectrum, r 1 and r 2 stand for the electrical field amplitude reflection coefficient of the front input and rear output cavity mirrors, respectively and φ k is the round-trip phase delay of the k-th order sideband. For a sideband of the q-th order, the power coupling efficiency is described by the first kind Bessel function J k−q (2β) of the (k−q)-th order, where β is the phase modulation index of the EOM. As there are two processes of phase modulation in a round-trip transmission inside of the linear cavity, the modulation index in Equation (2) is doubled as 2β. For an integrated ring cavity with a single section of phase modulation, the modulation index should be set as β [35,36]. It should be noticed that all the generated sidebands contribute to the k-th order sideband after a round-trip propagation. Therefore, the accumulation of all the coupled power from each sideband is required in Equation (2). At the same time, the power of the other sidebands is simultaneously changed together with the k-th order. To analyze the generated comb spectrum precisely, therefore, an iteration calculation of the inter-mode power coupling has to be applied. In fact, even the proposer M. Kourogi himself admitted that such an iteration calculation is too complicated [26].
In [30], an approximation calculation was made with less than ±15 adjacent sidebands, instead of all the sidebands ideally. The solution of a 1024 × 1024 sparse matrix equation was applied to describe the CEEOCG output spectrum. However, there are two problems with this solution. On one hand, there are obvious errors between the simulation and the experimental results for the higher order comb modes. On the other hand, the programming and processing of such an algorithm are still very time consuming.
To achieve a rapid characterization of the CEEOCG comb spectrum, M. Kourogi proposed another simplified exponential model in [26]. According to the transmission function of a standard Fabry-Pérot cavity [36][37][38], the electrical field of the transmission beam from a CEEOCG can be expressed as: where R stands for the equivalent power reflection coefficient of the CEEOCG cavity mirrors. It can be calculated as R = r 1 × r 2 . The item βsinω m t represents the phase modulation in a single-pass inside the CEEOCG. φ F stands for the residual phase delay in a single-pass inside the CEEOCG. After the detailed derivation with a series of mathematical approximations given in the Appendix A, the simplified model can be finally expressed as: Nowadays, the comb spectrum model in Equation (4) has become the most popular method for the spectrum characterization of CEEOCGs, including the microring resonatorbased CEEOCGs [35,36]. However, the approximations during the derivation of this simplified CEEOCG comb spectrum model cannot be fully met in reality. To compare with our proposed method, we would summarize the important approximation steps as follows.
Firstly, it requires that ω m t approaches zero to achieve the equivalent infinitesimal replacement of the sine item sinω m t as ω m t. For a commonly used phase modulator of 7~38 mm long and 9~40 GHz modulation frequency [28,39], however, the real value of ω m t is in the range of 1.87π to 2.28π. The residual phase will bring obvious error to the model.
More importantly, the derivation of Equation (4) requires the residual phase delay φ F to approach zero all the time. Only in this case, the exponent item exp[-j(φ F +βω m t)] can be equivalently infinitesimal, replaced as 1-jβω m t. This requirement severely limits the applicable field of the spectrum model. The real residual phase delay φ F is effected by the mismatch among the seed laser frequency, phase modulation frequency and the cavity resonance. As another factor that cannot be ignored, the intracavity material dispersion introduces an extra residual phase delay for the higher order comb modes as well.
Therefore, the existing simplified CEEOCG comb spectrum model can only be applied to certain conditions with limited accuracy. Many more cases with parameters of larger range variations cannot be simulated and characterized.

The Proposed Method for CEEOCG Spectrum Characterization
To derive an accurate and non-iterative method for the CEEOCG comb spectrum characterization, the electrical field of the laser beam was carefully analyzed during its propagation inside the CEEOCG cavity. According to Equation (3), the electrical field of the transmission beam from CEEOCG can be expressed as: where n stands for the n-th round-trip propagation inside the CEEOCG cavity, i.e., the count of the round-trip propagation. To simplify the exponential term in (5), the Jacobi-Anger identity was applied in the derivation [40]. Accordingly, the electrical field of the transmission beam from CEEOCG can be transformed as: where J k (x) stands for the first kind Bessel function of the k-th order. Assuming Equation (6) as a Fourier series, the transmitted electrical field of the k-th order sideband can be shown as: where β n = (2n + 1)β and φ Fn = (2n + 1)φ F stand for the phase modulation index and the residual phase delay of the n-th round-trip propagation inside of the CEEOCG cavity. It should be noticed that the transmitted electrical field intensity of the k-th order sideband is not related to the power coupling among sidebands in Equation (7). Instead, the irreversible increase of the propagation count n makes Equation (7) a non-iterative equation. To calculate the optical power intensity of the k-th order sideband, the electrical field of the k-th order sideband is multiplied by its conjugate as: The multiplication of the infinite plural terms in (8) can be classified into two categories. The first category contains all the multiplication of the same propagation count n. This multiplication process eliminates the exponent terms completely. The result is a summation of R 2n J k 2 (β n ) for n from zero to infinity. For the second category, the cross multiplication of different propagation count n should be calculated. To simplify the result, the products of the same difference of n are gathered. Accordingly, all the imaginary terms of the sine function are cancelled out. The result of cross multiplication for a certain propagation count n, named I tkn2 , can be expressed as: where m stands for the difference of the cross multiplication terms, which can be varied from zero to n. When we add up the results from both categories, the power intensity of the k-th order sideband can be expressed as: Inspired by the definition of the Matrix multiplication and Hadamard product of matrices, Equation (10) can be further simplified for computer simulation as follows: where (A Rn C Jk(0~n) ) stands for the Hadamard product of the A Rn and C Jk(0~n) vectors, which multiplies the two vectors element by element. The upper corner mark T is the symbol of vector transposition.
As the residual phase delay of a single-pass inside of the CEEOCG, φ F varies with the working condition of CEEOCG in real time. It consists of the mismatch phase delay φ α between the seed laser frequency and cavity resonance, the mismatch phase delay φ ∆f between the cavity resonance and phase modulation, and the phase delay φ D from the intracavity material dispersion. Ignoring the phase delay of 2nπ from even times of interface reflections, the residual phase delay of single-pass φ F can be expressed as: where ν FSR is the free spectral range (FSR) of the CEEOCG cavity, δν is the frequency difference between the seed laser and the adjacent cavity resonance, f m stands for the frequency of phase modulation, ω m = 2πf m . ∆f m is the mismatch between the phase modulation frequency and the cavity FSR, and GVD and L c stand for the group velocity dispersion and the length of the EOM crystal. According to Equation (17), the phase delay φ α , φ ∆f and φ D are all non-zero frequency dependent parameters in practical applications. The complete proposed method for CEEOCG spectrum characterization is composed of Equations (11)- (17). It should be noted that there is no more iterative calculation in this method. Once the parameters are determined, the comb spectrum of the CEEOCG can be calculated with a straight process. To characterize the CEEOCG comb spectrum with higher time efficiency, the vectors of A Rn and C Jk(0~2n) can be generated with the preset parameters and stored in advance. In this case, the calculation of (11)-(17) requires the summation of matrix multiplication and element-by-element multiplication only. All the processes can be easily and rapidly realized by the existing matrix computation software.
More importantly, it should be noted that there is no more mathematical approximation during the above derivation process. The proposed method for CEEOCG comb spectrum characterization is based on the accumulation of the Bessel function results with the phase modulation index β n of the round-trip propagation count n. Thus, the simulation accuracy of the proposed method is, in principle, determined by the maximum round-trip propagation count n max only. However, it should be noted that the truncation error of the floating-point numbers will influence the simulation accuracy as well, even for our non-iterative method. According to the ISO/IEC international standard 60559-2020 for floating-point arithmetic (i.e., IEEE standard 754-2019) [41], the minimum distance between two adjacent double-precision numbers is 2 −52 , i.e., approximately 2.220446 × 10 −16 . For the quadruple and even octuple precision floating point numbers, the truncation error is down to 2 −112 and 2 −237 , respectively. For the newly accumulated term ∆I tk of (10) with the increasing of the round-trip propagation count from n max to n max + 1, the attenuation property of the Bessel function ensures the decrease of ∆I tk to below this truncation error when n max is large enough. Hence, the convergence of the proposed method can be proved as well.
To analyze the influence of n max on the accuracy of the proposed method quantitatively, the spectrum of the center ±300 modes was simulated with n max = 100, 300, 1000 and 3000. The power reflection coefficient R and phase modulation index β were set to be 96% and 0.7 rad, respectively. When the residual phase delay φ F = 0, the simulated spectrum is shown in Figure 2a, where the simulated curves of n max = 300, 1000 and 3000 seem to be overlapped with each other. To ensure that n max is applicable for a more general case of non-zero φ F , the mismatch phase delay φ α was set to its lower limit of −0.7 rad. The mismatch frequency ∆f m and cavity FSR ν FSR were assumed to be 2 MHz and 9.2 GHz, respectively. Considering the GVD and length of LiNbO 3 EOM as 350.74 fs 2 /mm [42] and 10 mm, the spectrum simulation is shown in Figure 2b with the modulation frequency f m of 9.2 GHz. In this case, the curve of n max = 300 separates with the curves of n max = 1000 and 3000. Therefore, n max =1000 is large enough to simulate the center ±300 combs of a fixed CEEOCG in the working conditions above. To prevent redundant computations without improving the simulation accuracy, n max was set as 1000 for the following simulations of the CEEOCG spectrum with the power reflection coefficient R and phase modulation index β being 96% and 0.7 rad, respectively. down to 2 -112 and 2 -237 , respectively. For the newly accumulated term ΔItk of (10) with the increasing of the round-trip propagation count from nmax to nmax + 1, the attenuation property of the Bessel function ensures the decrease of ΔItk to below this truncation error when nmax is large enough. Hence, the convergence of the proposed method can be proved as well.
To analyze the influence of nmax on the accuracy of the proposed method quantitatively, the spectrum of the center ±300 modes was simulated with nmax = 100, 300, 1000 and 3000. The power reflection coefficient R and phase modulation index β were set to be 96% and 0.7 rad, respectively. When the residual phase delay ϕF = 0, the simulated spectrum is shown in Figure 2a, where the simulated curves of nmax = 300, 1000 and 3000 seem to be overlapped with each other. To ensure that nmax is applicable for a more general case of non-zero ϕF, the mismatch phase delay ϕα was set to its lower limit of -0.7 rad. The mismatch frequency Δfm and cavity FSR νFSR were assumed to be 2 MHz and 9.2 GHz, respectively. Considering the GVD and length of LiNbO3 EOM as 350.74 fs 2 /mm [42] and 10 mm, the spectrum simulation is shown in Figure 2b with the modulation frequency fm of 9.2 GHz. In this case, the curve of nmax = 300 separates with the curves of nmax = 1000 and 3000. Therefore, nmax =1000 is large enough to simulate the center ±300 combs of a fixed CEEOCG in the working conditions above. To prevent redundant computations without improving the simulation accuracy, nmax was set as 1000 for the following simulations of the CEEOCG spectrum with the power reflection coefficient R and phase modulation index β being 96% and 0.7 rad, respectively. To compare the existing approximation models based on the exponential function and the power coupling among ±15 adjacent sidebands with the proposed model, simulated spectrums and their deviations with the same parameters above are shown in Figure  3a,b, respectively. For the center ±300 modes, the overall comb spectrums of all three models in Figure 3a show an approximately linear power decay with the increasing of comb mode k. To compare the efficiency, the simulation time of the proposed model and the power coupling model were characterized by the matrix computation software as 1.2 s and 72.9 s, respectively. The detailed deviations are shown in Figure 3b with the proposed model as a reference. For the exponential approximation model, a linear error up to 1.6 dB can be observed for the higher order comb modes. Meanwhile, there is an extra 0.4 dB error for the ±1st comb modes. In contrast, the error of the power coupling model is nearly Figure 2. Simulated comb spectrums of the center ±300 combs with different maximum round-trip propagation count n max of 100, 300, 1000 and 3000. (a) Simulation with residual phase delay φ F = 0, (b) Simulation with φ F = −0.7 + 5.46 × 10 −3 k + 5.86 × 10 −6 k 2 rad, corresponding to the mismatch phase delay φ α = −0.7 rad, the mismatch frequency ∆f m = 2 MHz, the cavity FSR ν FSR = 9.2 GHz, the modulation frequency f m = 9.2 GHz, the GVD and length of LiNbO3 EOM being 350.74 fs 2 /mm and 10 mm, respectively.
To compare the existing approximation models based on the exponential function and the power coupling among ±15 adjacent sidebands with the proposed model, simulated spectrums and their deviations with the same parameters above are shown in Figure 3a,b, respectively. For the center ±300 modes, the overall comb spectrums of all three models in Figure 3a show an approximately linear power decay with the increasing of comb mode k. To compare the efficiency, the simulation time of the proposed model and the power coupling model were characterized by the matrix computation software as 1.2 s and 72.9 s, respectively. The detailed deviations are shown in Figure 3b with the proposed model as a reference. For the exponential approximation model, a linear error up to 1.6 dB can be observed for the higher order comb modes. Meanwhile, there is an extra 0.4 dB error for the ±1st comb modes. In contrast, the error of the power coupling model is nearly zero for the comb modes within the ±100-th order. With the further increasing of comb mode k, however, there is a rapid nonlinear raising of deviation up to 4 dB. We attribute this nonlinear error to the incomplete analysis of the power coupling model and the iterated accumulation of the floating-point truncation error. In summary, the exponential model shows a linear error but can only be applied when φ F = 0. The power coupling model can be applied with arbitrary φ F value but is not suitable for the analysis of CEEOCGs with abundant comb modes. zero for the comb modes within the ±100-th order. With the further increasing of comb mode k, however, there is a rapid nonlinear raising of deviation up to 4 dB. We attribute this nonlinear error to the incomplete analysis of the power coupling model and the iterated accumulation of the floating-point truncation error. In summary, the exponential model shows a linear error but can only be applied when ϕF = 0. The power coupling model can be applied with arbitrary ϕF value but is not suitable for the analysis of CEEOCGs with abundant comb modes.

Comprehensive Spectrum Characterization of CEEOCG with Different Parameters
Based on the proposed CEEOCG comb spectrum model above, the influence from the parameters of CEEOCG can be thoroughly characterized, e.g., the power reflection coefficient of cavity mirror R, phase modulation index β and phase delay ϕF. As the phase delay ϕF consists of the dispersion phase delay ϕD and the mismatch phase delay ϕα and ϕΔf, their independent and combined impact on the CEEOCG comb spectrum will be revealed in this section.

The Influence of Cavity Mirror Power Reflection Efficiency R and Phase Modulation Index β
As shown in Figure 4a, the power distribution of the generated comb spectrum was simulated with a power reflection efficiency R of 90%, 93%, 96%, 99% and 99.5%, respectively. In this simulation, the phase modulation index β and phase delay ϕF are assumed to be 0.7 rad and 0, respectively. According to the discussion above, the maximum propagation count nmax was set to be 1000 to simulate the center ±300 comb modes. With the increase of the power reflection efficiency R from 90% to 99.5%, the slope of power decay decreased from 0.653 to 0.006 dB per comb mode. Accordingly, the simulation curve became flatter. This simulation result fits the principle of the CEEOCG very well. The increasing of R introduces more power oscillation inside the cavity. Consequently, more power is distributed from the center to the higher order modes with the enhanced phase modulation.
When the generated comb spectrum was simulated with phase modulation index β of 0.3 rad, 0.5 rad, 0.7 rad, 1.2 rad and 2.5 rad, similar curves are shown in Figure 4b. In this simulation, the power reflection efficiency R and phase delay ϕF were assumed to be 96% and 0, respectively. The increasing of phase modulation index β from 0.3 rad to 2.5 rad leads to the decreasing of the power decay slope from 0.591 to 0.093 dB per comb mode. It corresponds to a flattening of the simulation curve and a broadening of the spectrum bandwidth. This phenomenon can be explained with an enhanced sideband generation capability during a single pass through the phase modulator. For the simulations with a much lower value of R or β, the sideband generation effect in a round-trip propagation through the phase modulator is much weaker. Thus, a larger propagation count nmax is required to simulate the accurate power of higher order comb modes.

Comprehensive Spectrum Characterization of CEEOCG with Different Parameters
Based on the proposed CEEOCG comb spectrum model above, the influence from the parameters of CEEOCG can be thoroughly characterized, e.g., the power reflection coefficient of cavity mirror R, phase modulation index β and phase delay φ F . As the phase delay φ F consists of the dispersion phase delay φ D and the mismatch phase delay φ α and φ ∆f , their independent and combined impact on the CEEOCG comb spectrum will be revealed in this section.

The Influence of Cavity Mirror Power Reflection Efficiency R and Phase Modulation Index β
As shown in Figure 4a, the power distribution of the generated comb spectrum was simulated with a power reflection efficiency R of 90%, 93%, 96%, 99% and 99.5%, respectively. In this simulation, the phase modulation index β and phase delay φ F are assumed to be 0.7 rad and 0, respectively. According to the discussion above, the maximum propagation count n max was set to be 1000 to simulate the center ±300 comb modes. With the increase of the power reflection efficiency R from 90% to 99.5%, the slope of power decay decreased from 0.653 to 0.006 dB per comb mode. Accordingly, the simulation curve became flatter. This simulation result fits the principle of the CEEOCG very well. The increasing of R introduces more power oscillation inside the cavity. Consequently, more power is distributed from the center to the higher order modes with the enhanced phase modulation.

The Independent Influence of the Mismatch Phase Delay ϕα, Mismatch Phase Delay ϕΔf and Dispersion Phase Delay ϕD
To reveal the influence of the mismatch phase delay ϕα, mismatch phase delay ϕΔf and dispersion phase delay ϕD independently, three simulations were performed with the same power reflection efficiency R of 96%, phase modulation index β of 0.7 rad and maximum propagation count nmax of 1000. Assuming the mismatch phase delay ϕΔf and dispersion phase delay ϕD to be zero, the CEEOCG comb spectrums with mismatch phase delay ϕα of 0, ±0.8β and ±β are shown in Figure 5a. In this situation, the symmetrical curves for different β of the same absolute value but opposite sign overlap with each other. With the increasing comb mode order, a linear power decay in the log scale can be observed. When the generated comb spectrum was simulated with phase modulation index β of 0.3 rad, 0.5 rad, 0.7 rad, 1.2 rad and 2.5 rad, similar curves are shown in Figure 4b. In this simulation, the power reflection efficiency R and phase delay φ F were assumed to be 96% and 0, respectively. The increasing of phase modulation index β from 0.3 rad to 2.5 rad leads to the decreasing of the power decay slope from 0.591 to 0.093 dB per comb mode. It corresponds to a flattening of the simulation curve and a broadening of the spectrum bandwidth. This phenomenon can be explained with an enhanced sideband generation capability during a single pass through the phase modulator. For the simulations with a much lower value of R or β, the sideband generation effect in a round-trip propagation through the phase modulator is much weaker. Thus, a larger propagation count n max is required to simulate the accurate power of higher order comb modes. To reveal the influence of the mismatch phase delay φ α , mismatch phase delay φ ∆f and dispersion phase delay φ D independently, three simulations were performed with the same power reflection efficiency R of 96%, phase modulation index β of 0.7 rad and maximum propagation count n max of 1000. Assuming the mismatch phase delay φ ∆f and dispersion phase delay φ D to be zero, the CEEOCG comb spectrums with mismatch phase delay φ α of 0, ±0.8β and ±β are shown in Figure 5a. In this situation, the symmetrical curves for different β of the same absolute value but opposite sign overlap with each other. With the increasing comb mode order, a linear power decay in the log scale can be observed. The linear slopes are 0.193, 0.421 and 0.806 dB per comb mode for φ α = 0, ±0.8β and ±β, respectively. The increasing of the mismatch phase delay φ α causes a concentration of the optical power to the center combs and there is a significate decrease in the CEEOCG comb bandwidth. This simulation result fits well with the previous reports in [33,34]. As shown in Figure 5b, the transmission rate of the CEEOCG varies with the mismatch phase delay φ α according to the theoretical analysis in [33,34]. The operation point of φ α = 0 enables a maximum power coupling from the incident laser mode to the higher order comb modes, which causes a minimum total transmission power. On the contrary, the operation point of φ α = ±β keeps most of the power in the lower order comb modes and achieves the maximum transmission power. in Figure 5d, respectively. The mode spacing is assumed to be 9.2 GHz. The dispersion phase delay ϕD is proportional to the square of the comb mode order k. Accordingly, the accelerated slope of power decay for the curves increases with the mismatch phase delay ϕD, as shown in Figure 5d. As the change of dispersion phase delay ϕD will not influence the total power inside of the CEEOCG cavity, the spectrum overlap of the comb modes within ±100 orders can be explained. From the analysis above, it is clear that the non-zero mismatch phase delay ϕα, ϕΔf and dispersion phase delay ϕD cause anarrowing of the comb spectrum. However, the corresponding comb spectrums are still symmetrical under the independent influence of all three phase delays. With the increasing of comb mode order k, only the slope of power decay for the mismatch phase delay ϕα is linear in the dB scale. For both the mismatch phase delay ϕΔf and dispersion phase delay ϕD, the power decay is accelerated. Meanwhile, only the mismatch phase delay ϕα will change the total transmission rate of the In Figure 5c, the influence of the mismatch phase delay φ ∆f is simulated by assuming that the mismatch phase delay φ α and dispersion phase delay φ D are zero. According to Equation (17), the mismatch phase delay φ ∆f is determined by the order of comb mode k, the mismatch frequency between cavity resonance and modulation frequency ∆f m and the FSR of CEEOCG cavity ν FSR . When the FSR of CEEOCG cavity ν FSR is assumed to be 9.2 GHz, the influence of the mismatch phase delay φ ∆f is simulated with a mismatch frequency ∆f m of 0, ±1 MHz and ±2 MHz. As the phase delay φ ∆f is proportional to the order of comb mode k, the decay of power is accelerated with the increasing of the comb mode order. The total power inside of the CEEOCG cavity does not vary with the mismatch phase delay φ ∆f . Therefore, the curves for the comb modes within ±50 orders are overlapped in Figure 5c. With a zero dispersion phase delay φ D , the symmetrical curves of the same absolute value of φ ∆f are consistent with each other as well.
In Figure 5d, the CEEOCG comb spectrums with different values of dispersion phase delay φ D are simulated by assuming the mismatch phase delay φ α and φ ∆f are zero. According to Equation (17), the dispersion phase delay φ D is mainly caused by the phase modulator. When the wavelength of the incident laser is 826.2 nm, the GVD value of the LiNbO 3 material can be found as 350.74 fs 2 /mm from [39]. Assuming the lengths of EOM crystal L c to be 0, 10 mm, 20 mm and 40 mm, the comb spectrum is simulated and shown as the black solid line, the blue dash line, the magenta dash-dot line and the red dot line in Figure 5d, respectively. The mode spacing is assumed to be 9.2 GHz. The dispersion phase delay φ D is proportional to the square of the comb mode order k. Accordingly, the accelerated slope of power decay for the curves increases with the mismatch phase delay φ D, as shown in Figure 5d. As the change of dispersion phase delay φ D will not influence the total power inside of the CEEOCG cavity, the spectrum overlap of the comb modes within ±100 orders can be explained.
From the analysis above, it is clear that the non-zero mismatch phase delay φ α , φ ∆f and dispersion phase delay φ D cause anarrowing of the comb spectrum. However, the corresponding comb spectrums are still symmetrical under the independent influence of all three phase delays. With the increasing of comb mode order k, only the slope of power decay for the mismatch phase delay φ α is linear in the dB scale. For both the mismatch phase delay φ ∆f and dispersion phase delay φ D , the power decay is accelerated. Meanwhile, only the mismatch phase delay φ α will change the total transmission rate of the CEEOCG. For lower order comb modes, the variations of the mismatch phase delay φ ∆f and dispersion phase delay φ D are not notable.

The Influence of the Mismatch Phase Delay φ α and φ ∆f with a Constant Dispersion Phase Delay φ D
The analysis above focuses on the independent influence of each phase delay. In a practical application of the CEEOCG, however, the dispersion phase delay φ D is usually a constant non-zero value. In contrast, the phase delays φ α and φ ∆f are always variable owing to the mismatch among the incident laser frequency, the cavity resonance frequency and the phase modulation frequency. Thus, the investigation of the influence of the mismatch phase delays φ α and φ ∆f with a constant value of dispersion phase delay φ D is of significance. The following spectrum characterizations are made with the same power reflection efficiency R of 96%, phase modulation index β of 0.7 rad and maximum propagation count n max of 1000.
Assuming the applied EOM is a 10 mm LiNbO 3 crystal with a GVD of 350.74 fs 2 /mm [39], the dispersion phase delay φ D can be calculated with Equation (17) as 5.86 × 10 −6 k 2 rad when the phase modulation frequency is 9.2 GHz. Assuming the mismatch phase delay φ ∆f is zero, the comb spectrums of the mismatch phase delay φ α = 0, ±0.8β and ±β were simulated and are shown in Figure 6a. If we compare it to Figure 5a, there are two significant differences. Firstly, the curves of the same absolute value of mismatch phase delay φ α are separated with different shapes. If we take the cases of φ α = ±β as examples, there is only a slight difference that can be observed for the center comb modes within the ±30 order. With the increase of the mode order k, the spectrum decays rapidly for φ α = β. For the phase delay φ α of −β, however, a spectrum broadening effect is shown for the combs from ±50 up to over ±300 orders. The newly generated comb modes exist as two low-energy wings of the spectrum. Thus, this phenomenon could be applied for the broadening of the CEEOCG spectrum with nearly no extra cost. The same situation could be found for the cases of φ α = ±0.8β as well. Secondly, the power decay of the curves is no more linear in the dB scale. The reason is that the introduction of the dispersion phase delay relates the total phase delay to the order of the comb mode. Of course, the symmetry of the curves in Figure 6a is kept the same as in Figure 5a. Besides, the concentration of optical power to lower order comb modes can be found for a larger absolute value of phase delay φ α as well.
simulated and are shown in Figure 6a. If we compare it to Figure 5a, there are two significant differences. Firstly, the curves of the same absolute value of mismatch phase delay ϕα are separated with different shapes. If we take the cases of ϕα = ±β as examples, there is only a slight difference that can be observed for the center comb modes within the ±30 order. With the increase of the mode order k, the spectrum decays rapidly for ϕα = β. For the phase delay ϕα of -β, however, a spectrum broadening effect is shown for the combs from ±50 up to over ±300 orders. The newly generated comb modes exist as two low-energy wings of the spectrum. Thus, this phenomenon could be applied for the broadening of the CEEOCG spectrum with nearly no extra cost. The same situation could be found for the cases of ϕα = ±0.8β as well. Secondly, the power decay of the curves is no more linear in the dB scale. The reason is that the introduction of the dispersion phase delay relates the total phase delay to the order of the comb mode. Of course, the symmetry of the curves in Figure 6a is kept the same as in Figure 5a. Besides, the concentration of optical power to lower order comb modes can be found for a larger absolute value of phase delay ϕα as well. As shown in Figure 6b, the influence of the mismatch phase delay ϕΔf with the same dispersion phase delay ϕD and zero phase delay ϕα was simulated. In these cases, the difference between Figures 5c and 6b can be summarized into two respects. Firstly, the symmetry of the curves themselves is broken by the mismatch phase delay ϕΔf with the help of the dispersion phase delay. The positive mismatch of Δfm leads to a left shifting of the simulation curve and vice versa, e.g., the blue dash line and the magenta triangle line in Figure 6b correspond to Δfm = 1 MHz and -1 MHz, respectively. The reason for this phenomenon is that the existence of the dispersion phase delay compensates for the positive mismatch of Δfm for the negative order of comb modes. At the same time, this effect makes the power of positive order modes decay faster. Of course, a constant dispersion phase delay can only compensate for the mismatch phase delay ϕΔf to a limited extent. Thus, the power decreasing of the negative order modes for the curve of Δfm = 2 MHz is no more linear as in the case of Δfm = 1 MHz, although the power decreasing rate of the negative order modes is always lower than the positive order modes. Secondly, the simulation curves of the same absolute value but opposite sign are symmetrical to each other. This phenomenon further proves that the change of curve shape from Figures 5c to Figure 6b is caused by a constant dispersion phase delay. At the same time, a reverse of the above phenomenon can be expected with an EOM crystal of negative GVD value.
A more comprehensive spectrum characterization should consider the variation of the mismatch phase delay ϕα and ϕΔf with a constant dispersion phase delay ϕD. By assuming the dispersion phase delay ϕD of 10 mm EOM being 5.86 × 10 −6 k 2 rad, the EOM modulation frequency and the FSR of CEEOCG cavity both being 9.2 GHz, the simulation results of the comb spectrum are shown in Figure 7. In Figure 7a, the mismatch frequency Δfm is set as the y-axis with scales of 0 MHz, ±1 MHz and ±2 MHz. The curves for different As shown in Figure 6b, the influence of the mismatch phase delay φ ∆f with the same dispersion phase delay φ D and zero phase delay φ α was simulated. In these cases, the difference between Figures 5c and 6b can be summarized into two respects. Firstly, the symmetry of the curves themselves is broken by the mismatch phase delay φ ∆f with the help of the dispersion phase delay. The positive mismatch of ∆f m leads to a left shifting of the simulation curve and vice versa, e.g., the blue dash line and the magenta triangle line in Figure 6b correspond to ∆f m = 1 MHz and −1 MHz, respectively. The reason for this phenomenon is that the existence of the dispersion phase delay compensates for the positive mismatch of ∆f m for the negative order of comb modes. At the same time, this effect makes the power of positive order modes decay faster. Of course, a constant dispersion phase delay can only compensate for the mismatch phase delay φ ∆f to a limited extent. Thus, the power decreasing of the negative order modes for the curve of ∆f m = 2 MHz is no more linear as in the case of ∆f m = 1 MHz, although the power decreasing rate of the negative order modes is always lower than the positive order modes. Secondly, the simulation curves of the same absolute value but opposite sign are symmetrical to each other. This phenomenon further proves that the change of curve shape from Figure 5c to Figure 6b is caused by a constant dispersion phase delay. At the same time, a reverse of the above phenomenon can be expected with an EOM crystal of negative GVD value.
A more comprehensive spectrum characterization should consider the variation of the mismatch phase delay φ α and φ ∆f with a constant dispersion phase delay φ D . By assuming the dispersion phase delay φ D of 10 mm EOM being 5.86 × 10 −6 k 2 rad, the EOM modulation frequency and the FSR of CEEOCG cavity both being 9.2 GHz, the simulation results of the comb spectrum are shown in Figure 7. In Figure 7a, the mismatch frequency ∆f m is set as the y-axis with scales of 0 MHz, ±1 MHz and ±2 MHz. The curves for different mismatch phase delays φ α of 0, ±0.5β and ±β are presented with a black solid line (φ α = 0), a red dash line (φ α = 0.5β), a blue dash-dot line (φ α = −0.5β), a magenta dash-dot-dot line (φ α = β) and a green short dash line (φ α = −β), respectively. Different from the analysis in Section 4.2, the mismatch phase delay φ α is changed from ±0.8β to ±0.5β for a better separation of the curves in Figure 7. In each slice of the joint analysis in Figure 7a, the distribution of each comb spectrum corresponds to the variation of φ α in sequence. For the slice of ∆f m = −2 MHz, the simulated curves are φ α = −β, −0.5β, 0, 0.5β, and β successively from left to right. The CEEOCG comb spectrum is not symmetrical for a certain non-zero mismatch frequency ∆f m . For the mismatch frequencies ∆f m of the same absolute value but opposite sign, however, a horizontal inversion of the simulated CEEOCG comb spectrums can be found. With the increase of the absolute value of ∆f m , the power decay of higher order comb modes is significantly enhanced. When the normalized power is less than −180 dB, some jitters can be found. We attribute this phenomenon to the floating-point truncation error during the simulation. For the unexpected extra spectrum in the modes over the ±250-th order, the cause of formation is still not clear. More work will be conducted on this in the future. Section 4.2, the mismatch phase delay ϕα is changed from ±0.8β to ±0.5β for a better separation of the curves in Figure 7. In each slice of the joint analysis in Figure 7a, the distribution of each comb spectrum corresponds to the variation of ϕα in sequence. For the slice of Δfm=−2 MHz, the simulated curves are ϕα = −β, −0.5β, 0, 0.5β, and β successively from left to right. The CEEOCG comb spectrum is not symmetrical for a certain non-zero mismatch frequency Δfm. For the mismatch frequencies Δfm of the same absolute value but opposite sign, however, a horizontal inversion of the simulated CEEOCG comb spectrums can be found. With the increase of the absolute value of Δfm, the power decay of higher order comb modes is significantly enhanced. When the normalized power is less than −180 dB, some jitters can be found. We attribute this phenomenon to the floating-point truncation error during the simulation. For the unexpected extra spectrum in the modes over the ±250-th order, the cause of formation is still not clear. More work will be conducted on this in the future. In Figure 7b, the mismatch phase delay ϕα is set as the y-axis with scales of 0, ±0.5β and ±β. The curves for different mismatch frequencies Δfm of 0 MHz, ±1 MHz and ±2 MHz are presented with a black solid line (Δfm = 0 MHz), a red dash line (Δfm = 1 MHz), a blue dash-dot line (Δfm = −1 MHz), a magenta dash-dot-dot line (Δfm = 2 MHz) and a green short dash line (Δfm = −2 MHz), respectively. Due to the spectrum-narrowing effect of the mismatch frequency Δfm, the distribution of CEEOCG spectrums in each slice of the same mismatch phase delay ϕα is not sequential with the increasing of Δfm. However, the horizontal mirror effect of the simulated spectrums from opposite mismatch frequencies Δfm is retained.
The simulation results above are not only a verification of the proposed CEEOCG comb spectrum characterization method. They can be further applied to identify the working condition of CEEOCGs. In our previous work [37], the pre-adjustment of a CEEOCG is very time consuming in order to match the incident laser frequency, the cavity resonance and the EOM modulation frequency at the same time. It was quite a hard process as the generated comb is very unstable. One has to be very experienced to be able to identify the working condition of the CEEOCG. Then, the further fine tuning can be implemented properly. With the proposed method and the simulation results, it will not be a In Figure 7b, the mismatch phase delay φ α is set as the y-axis with scales of 0, ±0.5β and ±β. The curves for different mismatch frequencies ∆f m of 0 MHz, ±1 MHz and ±2 MHz are presented with a black solid line (∆f m = 0 MHz), a red dash line (∆f m = 1 MHz), a blue dash-dot line (∆f m = −1 MHz), a magenta dash-dot-dot line (∆f m = 2 MHz) and a green short dash line (∆f m = −2 MHz), respectively. Due to the spectrum-narrowing effect of the mismatch frequency ∆f m , the distribution of CEEOCG spectrums in each slice of the same mismatch phase delay φ α is not sequential with the increasing of ∆f m . However, the horizontal mirror effect of the simulated spectrums from opposite mismatch frequencies ∆f m is retained.
The simulation results above are not only a verification of the proposed CEEOCG comb spectrum characterization method. They can be further applied to identify the working condition of CEEOCGs. In our previous work [37], the pre-adjustment of a CEEOCG is very time consuming in order to match the incident laser frequency, the cavity resonance and the EOM modulation frequency at the same time. It was quite a hard process as the generated comb is very unstable. One has to be very experienced to be able to identify the working condition of the CEEOCG. Then, the further fine tuning can be implemented properly. With the proposed method and the simulation results, it will not be a problem anymore. Based on the experimental setup established in Physikalisch-Technische Bundesanstalt [37], four CEEOCG comb spectrums were obtained and are shown in Figure 8. During the adjustment of the CEEOCG optical structure, these imperfect comb spectrums are present occasionally. With the dispersion phase delay φ D calculated as 5.86 × 10 −6 k 2 rad, the outer profile of the CEEOCG comb spectrums in Figure 8 is fitted with the simulation curves in Figure 7.  Figure 8d, respectively. By fine tuning the CEEOCG cavity length and the incident laser frequency according to the working condition identification above, the CEEOCG can be optimized to the ideal locking point in [37]. Therefore, the implementation of CEEOCG can be highly simplified with the help of the proposed comb spectrum characterization method. spectrums are present occasionally. With the dispersion phase delay ϕD calculated as 5.86 × 10 −6 k 2 rad, the outer profile of the CEEOCG comb spectrums in Figure 8 is fitted with the simulation curves in Figure 7. The positive and negative polarity of the mismatch phase delay ϕα and frequency Δfm can be first identified. Then, the very well-matched simulation curves can be obtained, with ϕα and Δfm being −0.98β and −110 MHz for Figure 8a Figure 8d, respectively. By fine tuning the CEEOCG cavity length and the incident laser frequency according to the working condition identification above, the CEEOCG can be optimized to the ideal locking point in [37]. Therefore, the implementation of CEEOCG can be highly simplified with the help of the proposed comb spectrum characterization method.

Conclusions
In summary, we proposed an accurate and comprehensive comb spectrum characterization method for CEEOCGs. The optical electrical field with respect to the count of the round-trip propagation inside of CEEOCGs was accumulated. The content of the residual phase delay calculation equation in the proposed method ensures its applicability to the arbitrary working conditions of CEEOCGs. The simplification of the proposed method was accomplished without any mathematical approximation by using the Jacobi-Anger identity and Euler's formula. With a maximum propagation count larger than 1000, the simulation error for the center ±300 comb modes comes from the truncation error of floating-point numbers only. Comparison results proved the error of the existing exponential approximation and the power coupling model to be linear and accelerated, increasing with the order of comb modes, respectively. Moreover, the proposed method can be efficiently computed with the Matrix multiplication and Hadamard product of matrices. By avoiding the iterative calculation of the power coupling among the generated comb modes, the simulation time can be reduced from 72.9 s down to 1.2 s.

Conclusions
In summary, we proposed an accurate and comprehensive comb spectrum characterization method for CEEOCGs. The optical electrical field with respect to the count of the round-trip propagation inside of CEEOCGs was accumulated. The content of the residual phase delay calculation equation in the proposed method ensures its applicability to the arbitrary working conditions of CEEOCGs. The simplification of the proposed method was accomplished without any mathematical approximation by using the Jacobi-Anger identity and Euler's formula. With a maximum propagation count larger than 1000, the simulation error for the center ±300 comb modes comes from the truncation error of floating-point numbers only. Comparison results proved the error of the existing exponential approximation and the power coupling model to be linear and accelerated, increasing with the order of comb modes, respectively. Moreover, the proposed method can be efficiently computed with the Matrix multiplication and Hadamard product of matrices. By avoiding the iterative calculation of the power coupling among the generated comb modes, the simulation time can be reduced from 72.9 s down to 1.2 s.
To reveal the influence of the parameters of a CEEOCG, a series of simulations based on the proposed method were conducted with the key parameters of CEEOCGs independently and jointly. The independent introduction of different non-zero mismatch and dispersion phase delays all led to an obvious narrowing of the comb spectrums. With the increase of the comb mode order, however, the speed of the symmetrical power decay is different for each type of residual phase delay. More realistic simulations were conducted by analyzing the two types of mismatch phase delays jointly with a constant dispersion phase delay. Hence, the exact power of each comb mode can be predicted accurately and efficiently. To the best of our knowledge, this is the first comprehensive characterization of all the key parameters of CEEOCGs.
The simulations and analyses above were not only a verification of the proposed CEEOCG comb spectrum characterization method. They can be further applied to identify the working condition of CEEOCGs. Four arbitrarily selected comb spectrums from an unstabilized CEEOCG were fitted with the proposed method. A high consistency can be found between the test data and the simulated fitting results, further proving the applicability and accuracy of the proposed method. With the help of fitting parameters for such a working condition identification, the CEEOCG can be further fine-adjusted or even optimally redesigned. Accordingly, the CEEOCGs with either a spatial linear cavity or an integrated ring cavity can service the applications in optical communications and waveform synthesis with a better performance. Therefore, we can find the electrical field intensity of the k-th order sideband from the first term of Equation (A2) as At the same time, the factor ω m t in the second term of (A2) can be obtained with the following property of Fourier series. ω m tg(ω m t) ↔ j dG k (k) dk . (A5) Consequently, the electrical field intensity of the k-th order sideband from the second term of Equation (A2) E tk2 can be expressed as As the power reflection coefficient R approaches 1 and the phase modulation index β is limited by the EOM to several rad, the intensity of the k-th order sideband in Equation (A6) is far smaller than the intensity in Equation (A4) and can be neglected. Considering the symmetrical power distribution of the positive and negative order of sidebands, the optical power intensity of the k-th order sideband can be finally approximated as where F is the finesse of the CEEOCG cavity. For R approaches 1, F ≈ πR/(1 − R).