Optimal Pulse-Shaping in Actively Q-Switched Ytterbium-Doped Fiber Lasers

In an actively Q-switched fiber laser (AQS-FL) a type of acousto-optic modulator (AOM) or (potentially) electro-optic modulator (EOM) controls the generation of output nanosecond wide pulses. An integrated Gaussian pulse shape is desirable in many applications such as material processing, microfabrication, ultrasound generation, gold photothermal therapy, etc. However, because of the system dynamics, generation of perfect Gaussian pulse shapes is not guaranteed in an AQS-FL, additionally designing the AQS-FL for a desired pulse peak and duration is an inverse problem which needs cumbersome trial-error efforts. We have developed a framework consisting of a rigorous FDM method plus a dedicated and innovative multi-objective genetic algorithm (GA) which assists the designer in achieving the desired Gaussian pulses within a reasonable time frame. The developed GA evolves the timing parameters of modulator plus the pump power and fiber length until the suitable goal is reached. To demonstrate the flexibility and design feasibility of our GA, three different single pulse and pulse train generation scenarios on a 7.5 m long Ytterbium-doped double clad fiber (YD-DCF) are examined to achieve the Gaussian,150 Wand 200 W peak power, 250 ns and 300 ns width pulses. To the best of our knowledge, it is the first implementation of an intelligent algorithm for optimizing the output pulse of an AQS-FL. It is worth noting that depending on the fiber host material and modulator specifications, much higher peak powers and different pulse durations are feasible, furthermore in case of utilizing the AOM, the pertaining limitations and feasibility are considered.


I. INTRODUCTION
A Gaussian pulse shape is easily defined by its peak power and full-width at half-maximum (FWHM) parameters, furthermore it causes an appropriate thermal profile in materials, as such the Gaussian profile pulses are desirable in different applications such as material processing, microfabrication, marking, medical equipment, LIDAR, etc. Actively Q-switched fiber laser (AQS-FL) is a major candidate for the above applications and is able to generate single or train of nanosecond wide pulses using either an acoustooptic or electro-optic modulator (AOM or OEM). However, a regular pulse shape or desired repetition regime under the applied switching frequency is not always guaranteed in AQS-FLs and the system may tend to produce non-Gaussian or multi-peak pulses or would be attracted to an unwanted repetition mode (attractor) [1]. In this work we introduce a method which finds the best possible combination of The associate editor coordinating the review of this manuscript and approving it for publication was Huaqing Li . parameters to avoid the time-consuming process of trial-error and minimize the AQS-FL design time.
AOMs with less than 10 ns transition times (e.g. TeO 2 devices) are available today [2], in addition the availability, damage threshold and control units of EOMs are under constant improvement, for example the damage threshold in some Rubidium Titanyl Phosphate (RTP) EOMs reaches to 1 GW/cm 2 [3]. In case of using the AOM, there are two inherent limitations, firstly the AOM rise time (t r ) and fall time (t f ) are limited by the sound speed and beam diameter, secondly, the rise and fall times should be assumed as equal parameters. We have taken into account these limitations when an AOM is applied in our algorithm. The on time (t on ) of AOM can be controlled electrically while t r and t f can be adjusted optically just in some extent, hence we aim only at selecting an appropriate and available AOM. Since AOM application in AQS-FL is more common, in this paper we have utilized an EOM only for conceptual (not practical) representation of our method in one case.
The methods of Q-switching operation can be classified into the active and passive switching. In the passive Q-switching (PQS), a saturable absorber (SA) like the Samarium-doped fiber or the dye SA performs the function of Q-switching, whereas in AQS this function is carried out by a type of modulator such as AOM or EOM. PQS takes advantage from the simple design, light weight, low cost and more consistency with the classical Q-switching theory, while the AQS is more advantageous from the viewpoint of precision and controllability. However, AQS includes a few more challenges such as the generation of multi-peak or separated output pulses with unwanted and hardly predictable overall shapes. This behavior of AQS has been referred frequently as the multi-peak phenomenon (MPP) in the literature [4]- [7]. Wang et al. in [4] explain the MPP using perturbation concept. Upadhyaya et al. in [7] introduce the Brillouin scattering as a major contributing factor to the MPP generation. Barmenkov et al. experimentally explored the repetitive regimes of AQS in the context of chaos theory [1]. In addition, other numerous literatures including [5], attempted to explain the physical origins of MPP by tracing the dynamics of system or suggested some modifications to the laser cavity using extra components [8]. Amongst all of the methods, modification of modulator timing parameters such as the rise time (t r ), on time (t on ) and fall time (t f ) ( Fig. 1) through trialerror, had been the most tried and feasible solution. The aim of designing an AQS-FL is to generate the output pulses having four desired properties including the shape, peak power (or energy), FWHM and repetition regime. Trial-error and analytical methods can be considered as two design solutions; however, the chaotic, time varying and nonlinear behavior [1], [5], [9] of AQS-FL systems, renders both of these methods unsuitable. On the other hand, since this is an inverse problem, an artificial intelligence (AI) method such as genetic algorithm (GA) seems promising [10], [11]. As such, we have developed and implemented a dedicated multi-objective GA to find the best (globally optimized) combination of parameters which is required to produce the MPP-free Gaussian pulse shapes with the desired peak power and duration. Taking benefit from intelligent methods such as GA, swarm or animal behavior algorithms has been carried out already on some types of laser systems [10], [12]. However, to the best of our knowledge, the present work is the first implementation of such algorithms on AQS fiber laser systems. Our developed GA code features some innovations to fulfill the consistency with fiber laser simulation methods, while its convergence speed is also improved using both parallel processing and AI techniques.
Here, a linewidth of 2 nm is considered, additionally the demonstrated numerical simulations are chosen to be under 200 W peak power in order to avoid the nonlinear effects such as stimulated Brillouin and stimulated Raman scattering (SBS and SRS) in worst case scenarios where unrealistically high nonlinear gains for the host glass are assumed. In practical cases however, the thresholds for this nonlinear effect are comparatively higher and there are several mitigation techniques which enable the designer to escalate the output peak power dramatically [13], [14]. The addition of different ceramic materials to the fiber core (e.g. in aluminosilicate fibers) also provides a promising method to decrease the nonlinear gain coefficients [15].
The paper is organized as follows. In section II, the theoretical background of AQS-FLs is summarized. Section III is dedicated to the numerical modeling. In section IV, our proposed GA method has been described in details. Section V consists of simulation results and discussion, and finally, section VI concludes and summarizes the paper.

II. THEORETICAL BACKGROUND
A typical configuration of Ytterbium-doped AQS-FL is shown in Fig. 2, where the forward pumping is performed using a single or an array of laser diodes and a modulator (AOM or EOM) performs the function of Q-switching element. In the continuous wave (CW) regime, noise helps to initiate the laser oscillations and amplified spontaneous emission (ASE) can be considered as a type of noise which is involved in the above process. ASE thereafter loses its major role to main lasing process and would be regarded as a parasitic phenomenon degrading the gain and efficiency of the system. Whereas in pulsed regimes such as Q-switching, the system is always engaged in the inception process of lasing which is the dominance scope of ASE. For CW simulation of fiber lasers, one can even ignore the effect of ASE and deploy just an approximate model, while in Q-switching, ASE would be the major theme of our modeling and simulations.
The gain coefficient of an active fiber cavity can be expressed as [16] g(z, λ k ) = k (λ k )[σ e (λ k )N 2 (z) − σ a (λ k )N 1 (z)] (1) where λ k is the k-th wavelength component, k is the optical overlap factor, N 1 and N 2 are the populations of the first and second energy levels, respectively, and σ e (σ a ) is the emission (absorption) cross-section of rare earth (RE) ions (Yb 3+ ). For an active fiber with the length of L, the gain of medium is defined as [16] where α is the fiber attenuation coefficient at the signal wavelength and g is defined in (1). Quantum defect (QD), the difference between pump and emission energies, not only defines the amount of heat dissipated in the cavity and hence, efficiency of the system, but also gives a clue to ASE spectral range. QD for Ytterbium-doped fibers is around 9% [17], which is considerably lower than the same parameter for other types of RE dopants like Neodymium (Nd) or Erbium (Er). As such, the pump and lasing wavelengths are comparatively near in Ytterbium-doped system. Pumping is done in 900 nm range (typically 915 and 976 nm) and ASE emission begins around 1020 -1030 nm and extends to around 1100 nm. In this context, modelling of gain and ASE spectra often is done for the same range [18].
Upon Q-switching, different ASE spectral components begin to compete, dedicate and consume their share of population inversion and available gain. A few of ASE components which obtain most of the share, grow much more rapidly than others, and most of the energy stored in the gain medium would be dedicated to them [18].
Relaxation oscillations (RO) in a laser system can be simply described as the transient oscillations after pump switch-on and before reaching the steady state. RO should be taken into account if a very rapid Q-switching is carried out. However, in AQS, the modulator transmission ramps up steadily and gradually, so the cavity deals with a type of smooth switching effect [4], [5]. Study of RO effect on the AQS output pulse shape would be possible through switching on the pump power just in the time of modulator ramp up under different rise time (t r ) conditions. Furthermore, a realistic and precise version of the absorption and emission cross-sections would be necessary to achieve correct and trustable results. Approximation of the cross-sections using Gaussian expansion (inhomogeneous broadening of host glass) is possible. The above expansion can be performed around transition wavelengths corresponding to Ytterbium upper and lower level manifolds. The group of resulting Gaussian profiles then constitute the corresponding cross-section curve [19], [20]. Additionally, one can also take into consideration the homogenous broadening of host material and obtain the convolution of Gaussian and Lorentzian distributions as the so-called Voigt profile [21].
However, utilizing the curves offered by the original equipment manufacturer (OEM) experimental data would yield the most accurate and realistic results. This data is usually available as a matrix consisting discretized samples of wavelengths, absorption and emission cross-sections, which covers the entire range of relevant spectrum.

III. NUMERICAL MODEL AND SIMULATION
The basis of AQS simulation is the following set of coupled partial differential equations (PDE) including population rate equations (3) -(4) and propagation equations (5) -(6). The rate equations explain the dynamics of upper and lower level populations, and the propagation equations express the evolution of pump and signal powers in both the forward and backward directions [5], [9], [22]- [25]. The pump light is assumed to be distributed homogenously across the inner cladding.
N 0 is the doping density of active fiber, N 1 and N 2 are the populations of upper and lower levels, P pf (P pb ) is the forward (backward) pump power and P kf (P kb ) is the forward (backward) signal power. The effect of ASE spectral shape is realized by sweeping the ASE wavelength (λ k ) where the index k represents the selected so-called ASE channel and K is the total number of channels. The last term in the right-hand side of equation (6) denotes the effect of spontaneous emission on the mode, in this numerical model a non-polarized state is assumed and the coefficient of 2 is due to this assumption [5]. Ytterbium-doped gain media have a quasi-three level system in which the upper level lifetime is around one millisecond and assuming a simple two-level system is adequate in almost all cases for modeling the Q-Switched YDFLs, where the energy levels can be divided into lower and upper state levels with N 1 and N 2 populations [21], [26].
In our implemented model, Rayleigh back-scattering (RBS) [27] is taken into account through the Sα RS P kb(f ) term in (6). If the difference between refractive indices of the core and inner cladding is small, for example if the numerical aperture (NA) is around 0.06, then it would be possible to consider the same velocities for the pump and signals (v p = v k ).
Descriptions, common values or ranges of other parameters in the above equations can be found in Table 1 [5], [9], [18], [22]- [25]. We have used different values of every parameter according to the ranges in Table 1 for performing numerous simulations runs and making sure about the validity of the developed numerical model. However, to utilize this model in the described GA method (section IV), the simulation setup details and selected values of parameters are explicitly specified (section V). The above set of coupled PDEs can be solved using different finite difference method (FDM) schemes; however, the variation of N 1 and N 2 along the considerable length of the fiber cavity should be addressed. Hence a traveling wave method (TWM) [22]- [25], [28]- [30] plus a stable forward time centered space (FTCS) combinational FDM scheme is applied through the programming code to calculate the population density, pump and signal powers in every 10 cm of the fiber length. In order to achieve a valid result, the simultaneous effect of different ASE channels (competition of spectral components) should be calculated at every specific time and position step. This process is accomplished via solving a number (K ) of the above PDE sets simultaneously or through a self-consistent iteration loop; otherwise the results would not be trustable. In this manner, by implementing the temporal, spatial and spectral discretized steps ( t, z and λ), a detailed spatial and spectro-temporal picture of the system conditions would be realized.
Boundary conditions (BCs) are imposed for the pump and laser powers at the beginning and end of the fiber cavity [9].
The parameters used in (7) - (9), their description and common values or ranges [5], [9], [18], [22]- [25], are given in Table 2, additionally these BCs are explicitly specified for the simulation model setup in section V. All of these parameters except the pump power are unitless. The R OC and R HR parameters can be assumed constant or wavelength dependent in the numerical model. In the case of using the wavelength selective mirrors or FBGs, one or a limited number of ASE components (λ k ) will be selected for lasing. In the case of single wavelength selection by FBGs, this wavelength would be denoted by λ s . Finally, according to Fig. 2, the output pulse power at the OC end of the fiber (z = L) can be expressed as [5] P out (L, t) =

IV. THE DEDICATED GENETIC ALGORITHM
We have developed and implemented a dedicated GA to approach a multi-objective goal. Here, the rise time (t r ), on time (t on ) and fall time (t f ) of the modulator transmission function plus the pump power and fiber length are regarded as chromosomes. The versatility of the proposed GA also makes it possible to extend the algorithm to different parameters such as doping density or fiber core diameter in order to optimize the system components. In case of utilizing an AOM, t r = t f and we aim at selecting the proper AOM rather than adjusting its timing parameters.
Each set of the above parameters would be an individual and a number of individuals constitute a generation. The first generation is randomly generated within a reasonable and broad range for each chromosome. A Gaussian goal function with a desired peak power and full-width at half-maximum (FWHM) is defined initially by the user. Next, the generated output pulse by every individual in the generation is compared with the above Gaussian and a fitness value is assigned to its fitness chromosome. Then the individuals in each generation would be sorted according to their fitness, inspected and classified into direct, crossover, mutation and deletion classes through marking their privilege chromosome.
Individuals with direct privilege are then transferred directly to the next generation, while the ones within crossover class go through a random crossover process along with already selected individuals. Furthermore, individuals marked with mutation privilege pass a mutation process before their crossover. Fig. 4 describes the flow chart of the GA method. The mutation process here is inspired by animal's behavioral models such as the one described in [31]. The mutation rate is adjusted dynamically to ensure the algorithm wouldn't be trapped in one of the local optima. Finally, the individuals within the deletion class undergo the extinction process in which they will be completely randomized again and can't find a path to the next generation. The described process is repeated in a loop until the fitness value is converged to the desired result. Each generation is created within one of the GA loops. Once the first round of the GA is converged, the resulted rise and fall times can be compared to the commercially available AOMs and the AOM with nearest specifications would be selected. Next, the rise and fall time chromosomes are set to a fixed value from the chosen AOM data and a second round of GA is started, in this second run, the fiber length parameter can be enabled as a chromosome or the GA can converge using only the evolution of on time and pump power. The second round of the GA would be much shorter since plenty of AOMs are available today, thus the result of the first GA round would be near to one of them just with few percentages of difference. A same procedure is performed for one of the simulation cases of section V to demonstrate the feasibility of the method.
In comparison to GA, some advantages of particle swarm algorithms (PSO) are including the simplicity, ease of development and higher convergence speed while GA offers better escaping opportunity from local optima and is easier to understand [32]. In our developed algorithm we have implemented a customized process for the mutation step; here the different chromosomes have not equal mutation rate (probability), for example the modulator fall time features the lowest mutation probability because typically the fall time of the modulator is not so effective in changing the pulse shape or pulse specifications, the mutation rate for each chromosome can be defined separately which gives a higher robustness to the GA. In addition, the developed GA features parallel processing which partially compensates the lower convergence speed in comparison to PSO algorithms.

V. SIMULATION MODEL SETUP, RESULTS AND DISCUSSION
The goal of the present work is to reach a desired peak power, FWHM and Gaussian pulse shape either in single pulse or pulse train generation, in one case we have also demonstrated the minimization of the fiber length to reduce the system cost. Executing the GA with 300 generations on a 4 cores Core i5 PC with 8 GB of RAM took approximately 8 hours. The simulation run time is inversely proportional to the number of CPU physical cores and directly proportional to the number of generations.
Based on the schematic of Fig. 2, a 7.5 m long silica YD-DCF is utilized as the laser cavity, however in some cases the fiber length would be variable and we intend to change or minimize it. The chosen spatial step size ( z) is 5 cm and it varies in case of fiber length optimization. Doping concentration is assumed to be 5.5 × 10 25 ions/m 3 and hence no ion-ion interaction is considered in the simulation [33]. Core/inner cladding diameters would be 10/125 µm and NA equal to 0.06 is assumed, thus the system operates in single transverse mode. Typical timing parameters in the range of tens and hundreds of nanoseconds have been set for the modulator operation and the pump power range of 5 -30 W has been considered. Sinusoidal ramp-up and ramp-down similar to Fig. 1 are set for both edges of transmission function. Switching frequency of 40 kHz has been set.
The pump and signal overlap factors ( p and k ) and loss coefficients (α p and α k ) are calculated based on the geometry of cylindrical DCF and their wavelength dependency is considered in the model [26]. Two FBGs according to Fig. 2, are adjusted in the numerical model to maintain the lasing spectrum around 1080 or 1064 nm which are the most well-known and essential wavelengths in different industrial, medical, scientific and other applications. For this purpose, R HR is assigned to be a constant value equal to 0.995 and R OC is set to 0.04 for the lasing wavelength. Other BCs were assigned the following values, η sp = 0.98, η p = 0.85, η k = 0.7. In order to keep a balance between the accuracy and total simulation time, a reasonable number of ASE channels (K = 36) are modeled. The modeling is based on the absorption and emission cross-sections of Ytterbium-doped fibers [33]. Fig. 5 illustrates three situations of the output pulse with arbitrary picking of t r , t f , t on and P p parameters. The curves for cases ii and iii are scaled down to create a proper visualization in this figure. The pulse shape for the cases (i) and (ii) of Fig. 5 is not regular, so it is not possible to aim at the parameters such as pulse FWHM. Furthermore, it would not be possible to define the pulse peak and duration using the classical solid state laser theory, which is just partially applicable in fiber lasers. In the case (iii) of Fig. 5 the pulse shape is regular; in this case, t r is dramatically increased in order to generate an MPP-free output pulse. However, designing an AQS system to create a predefined pulse peak power and desired FWHM isn't easily realizable.
Here, two Gaussian pulses with 150 W and 200 W peak powers and the FWHM of 250 and 300 ns are set as the goal function of our GA with 50 individuals per generation. After around 300 generations, the pulse specifications reach to the vicinity of the defined goal with ∼5% fitness error and after 1000 of generations the goal would be achieved with 1% -2% fitness error. The convergence of the GA gradually slows down with the increase of generation number. Increasing the mutation and randomization rates provides the escaping opportunity from the local optima; however, this should be limited to avoid a very slow convergence at the same time. Fig. 6 demonstrates the optimized output pulses for the two above mentioned cases after 600 generations. Fig. 6(c) shows the result of a second GA round which is performed for the case of Fig. 6(b), in Fig. 6(b) the AOM rise and fall times are equal to 335.88 ns which is near to the same parameters for an available AOM with 350 ns rise and fall times, so in second round, t r and t f were set to 350 ns and the GA started with (without) the fiber length parameter as one chromosome. The resulted values are presented in the following and the second GA round realizes the feasibility of our method after just tens of generations. The same goal can be reached with t r = t f = 350 ns, t on = 1737.3 ns, L = 7.09 m and P p = 21.71 W in case of enabling the fiber length chromosome for the second GA round.
A 3D fitness map resulting from running the GA with 300 generations is shown in Fig. 7, where mutation and randomization are visible through some dips inside this figure. Increasing the number of generations leads to more and wider dip areas, and correspondingly, the fitness value of the elite  individuals would be improved, which is an indication of the implemented elitism in our algorithm. Fig. 8 illustrates the generation of a Gaussian pulse with a peak power of 200 W and duration (FWHM) of 250 ns utilizing the above mentioned setup and different FBG central wavelength (1064 nm). Comparing with the two previous simulations at 1080 nm (Fig. 6), a small modification has been done in the model (R OC parameter in Table 2) to change the lasing spectrum; furthermore, the pump power is increased around 20% by the GA. The inset of Fig. 8 shows that the spectrum is shifted to 1064 nm lasing wavelength. In addition, the normalized population curve of the excited level is inserted into this figure which displays less than 6% of excitation. Fig. 9 illustrates the result of our GA optimization for the repetition regime of the same AQS-FL which is described  in simulation setup. The switching frequency is 40 kHz and sampling are done on two consecutive pulses; the system reaches to the desired state after generating a few pulses. Similar to Fig. 6(b) the peak power and FWHM are set to 200 W and 300 ns. The corresponding global optimum parameters are found to be t r = t f = 388.4 ns, t on = 2398.3 and Pp = 6.36 W using AOM. The rise and fall times can be set to the nearest available value (e.g. 350 ns) similar to what is performed for Fig. 6(c) and a second GA round will result in a practically feasible goal. Fig. 10 displays the result of fiber length optimization while we aim at the pulse peak power and FWHM of 200 W and 300 ns. In the previous optimization cases, we had considered a fiber length of 7.5 m while it's possible to reduce the fiber length in some extent. The results of the GA show that the fiber length can be reduced to 5.91 m which in turn reduces the cost of the optimized system. Fig. 11 shows the variations of t r , t on , t f and P p during the GA process for the simulation conditions corresponding to Fig. 6(a) and utilizing an EOM, it's supposed that the timing parameters t r and t f can be varied freely and independently. It is worth pointing out that having the values greater than 20 ns for these two parameters is not practically common in EOMs, thus this figure illustrates a conceptual representation of the GA evolution. The individual with the highest fitness value in each generation has been selected for this purpose.  As can be seen in the chart, the pump power only accepts some small variations during the GA after a few numbers of generations. In fact, the pump power can be altered in a wide range without changing the situation of MPP correction. However, this measure alters the peak power and FWHM of the output pulse. Therefore one can utilize this scalability of the system in its optimized state to change the pulse peak power and duration without losing the Gaussian shape.

VI. CONCLUSION
In summary, we have developed a dedicated finite deference method model based on the distributed traveling wave model for actively Q-switched fiber lasers (AQS-FLs). Furthermore, we have developed a dedicated genetic algorithm (GA) from scratch to optimize the output pulse and repetition characteristics of the AQS-FL. For this purpose, four optimization objectives including the MPP-free Gaussian pulse shape, desired peak power, pulse duration (FWHM) and repetition mode have been taken into account to find the best combination of the input parameters i.e., the rise time, on time, and fall time of the modulator as well as the pump power. Optimization is done for two versatile wavelengths of 1064 and 1080 nm. We evaluated the performance of the developed GA for optimizing an AQS-FL with a 7.5 m long YD-DCF and the numerical aperture (NA) of 0.06. The shape of the output pulse has been modified to have a Gaussian profile with the desired peak powers of 150 W and 200 W and the pulse durations of 250 ns and 300 ns in either single pulse or pulse train generation to demonstrate some samples of the feasibility and versatility of the proposed GA. Our simulation results demonstrate that the developed GA can successfully optimize the output pulse regarding its shape, peak power, duration and repetition mode.
The proposed methods of this paper pave the way for intelligent pulse shaping in AQS-FLs and are extendable to much higher peak powers depending on the host glass specifications. The method is even extendable to other gain media and lasing wavelengths in case of accounting for the variations in energy level structure.