Supercontinuum optimization for dual-soliton based light sources using genetic algorithms in a Grid platform

We present a numerical strategy to design fiber based dual pulse light sources exhibiting two predefined spectral peaks in the anomalous group velocity dispersion regime. The frequency conversion is based on the soliton fission and soliton self-frequency shift occurring during supercontinuum generation. The optimization process is carried out by a genetic algorithm that provides the optimum input pulse parameters: wavelength, temporal width and peak power. This algorithm is implemented in a Grid platform in order to take advantage of distributed computing. These results are useful for optical coherence tomography applications where bell-shaped pulses located in the second near-infrared window are needed.


Introduction
The Soliton self-frequency shift (SSFS) [1,2] plays a central role in many effects taking place during supercontinuum (SC) generation in optical fibers (see Refs. [3,4] for a review on the topic). To mention only a few examples, light trapping [5], multi-peak soliton states [6][7][8], emission of Airy waves [9], intense dark-soliton SC [10] or broad and intense blue shifted polychromatic dispersive waves [11,12], would not be possible (or strong enough) without the SSFS. One of the most notorious feature of the Raman effect in SC generation with femtosecond pulses corresponds to the Raman soliton carrying the lowest frequency. Its large frequency shift from the laser pulse has motivated infra-red (IR) Raman soliton sources [13][14][15] and their optimization [16][17][18][19][20].
In the previous work of Ref. [12], we have shown that fs-pulses traveling in a dispersion engineered single mode fiber (SMF) can generate several pre-defined spectral peaks in the normal group velocity dispersion (GVD) region. These peaks correspond to the narrow band Cherenkov radiation emitted by the bright solitons [21] undergoing Raman red-shift and recoil. Such spectra were important for applications in optical coherence tomography (OCT) with wavelengths in the near infra-red (NIR) window [22][23][24]: λ 1 µm.
In the present work, we are interested in OCT applications enabled by fiber based illumination [25,26] in the second near IR window (NIR II) [27][28][29] where λ ∈ [1, 1.4] µm, typically in the anomalous GVD regime of highly nonlinear photonic crystal fibers (PCFs). Because of the typical dispersion landscape, sources in the NIR II window may be based on the bright optical solitons arising during SC generation through the intricate soliton fission effect [30,31]. With this picture in mind, we used an efficient and general computational optimization method based on genetic algorithms (GA) [32] that is capable of finding output spectra [see Fig. 1] exhibiting two peaks centered at pre-defined wavelengths. This type of multi-peak spectral illumination is very often required for OCT applications in the NIR II [29] and the simultaneous presence of the two operating spectral components is a prerequisite for real time imaging [33,34]. The two peaks we obtain are the first and second Raman solitons presenting a clean bell-shaped spectral profile, essential for OCT [35], with widths providing a decent longitudinal resolution l c ≈ 10 µm [36]. To demonstrate the usefulness of this method, we consider in this work the two spectral channels separated by 100 nm (see, e.g., Ref. [29]). The optimization method finds the optimal input pulse parameters, namely central wavelength, λ 0 , temporal width, T 0 , and peak power, P 0 , yielding the desired spectra. The obtained peak powers are of up to 90 mW for each spectral band, satisfying the needs for OCT imaging applications [35]. We find by this method the possibility to tune the wavelength of the target spectral channels, what represents an important feature of this strategy since the greatest potential of the given PCF is exploited, specially in situations where limited choice of PCF designs is available.
It is worth mentioning that the use of GAs in optics has indeed proofed useful previously in solving satisfactorily the inverse optimization problem of the one we present here, i.e., the design of PCFs to control SC dynamics, in a wide range of situations [37][38][39][40][41]. The use of GAs generally requires a large amount of simulations. For this reason we used a distributed computing (Grid) platform to reduce the time required to find the optimal solutions. The advantage of this infrastructure is that it enables the use of the same code in a platform of scalable resources which are adapted according to the needs of the particular problem.

Supercontinuum modeling and the Genetic Algorythm
We simulate the nonlinear propagation of the complex electric field envelope, A, along the fiber axis, z, by integrating numerically (with split-step fourier method) the generalized nonlinear Schrödinger equation (GNLSE) [42], where the β q 's account for the linear fiber dispersion and is the nonlinearity parameter [43] which has been computed with a FEM solver (comsol) by integrating the electromagnetic components of the modal field at λ = 800 nm along the transverse fiber cross section. ε = 2.09 and n 2 = 2.6 × 10 −20 m 2 /W are the relative permitivitty and nonlinear index of silica glass, respectively. The nonlinear response of the glass is [44], where δ (t) is the Dirac delta function and the Raman (delayed) contribution is weighted with f R = 0.18 and described by the damped oscillator where τ 1 = 12.2 fs, τ 2 = 32 fs, and Θ(t) is the Heaviside step function. The input pulses used in our simulations are of the form √ P 0 sech(t/T 0 ), where (for a fixed λ 0 ) the input peak power controls the soliton order N ≡ T 0 γP 0 /|β 2 |.
We consider a length L = 25 cm of the NL-2.4-800 PCF, exhibiting the lowest zero GVD at 798 nm [see details in Fig. 1], convenient for broad band SC generation with Ti:Sapphire lasers. The frequency conversion performance is investigated within the attainable ranges of input pulse parameters: λ ∈ [750, 850] nm, P ∈ [5,15] kW, and T ∈ [30,150] fs. For each simulation along the PCF, the GA generates an individual with the genome | g ≡ [g 1 , g 2 , g 3 ] T = [T 0 , λ 0 , P 0 ] T and evaluates how suitable that individual is from the simulation output through the fitness function (to be minimized) defined as where 2∆ω is the chosen spectral channel widths and ω c j = 2πc/λ c j the central frequency.
Note that the definition of φ as a product tends to favor output spectra in the form ψ 1 ≈ ψ 2 600 800 1000 1200 1400 1600 1800  Fig. 2(a). After an initial set of p th randomly (uniformly distributed) generated individuals (stage 1), the genetic operators (GOs) crossover,X , and mutation,M , are responsible for generating the new offsprings, which are added to the population until the maximum size p = p max > p th is reached (stage 2). During stage 3 the population size is kept constant, p = p max , and a replace the worst strategy is used, i.e., the offspring is added to the population if φ o f f spring < φ max or disregarded otherwise. This steady state GA keeps the Grid constantly computing new individuals in parallel and fully exploits the processing power of the Grid. We briefly describe the GOs below.
Cross-over,X , is the first GO applied to the current population and generates two childs | g 1,2 c by combining two randomly chosen parents, | g 1,2 p , i.e.,X [| g 1 0 the empty set). We have used simulated binary crossover (SBX) [45] and the [12 × 12] operator where0 ≡ 0×Î (beingÎ the 2×2 identity matrix). The crossover activators, x k , set a probability for cross over of 95% per gene. The stochastic variables in this case,σ k [see Fig. 2(b)], are chosen from a uniform random number u k ∈ [0, 1] according to, The polynomial mutation,M , [46] suitable for real coded problems (continuous valued variables), is applied afterX and generates new genes asM :| g →| g ′ , where 2∆ k is the interval size of each variable (∆ τ = 60 fs, ∆ λ = 50 nm, ∆P = 5 kW). Hence, in average only one gene is mutated per individual when mutation is applied: with n = 20, peaked around ζ = 0 and clearly different from the random generation. For each gene the stochastic variableζ is chosen from a random u ∈ [0, 1]: Statistically,M provides diversity to the population andX explores the parameter space in the vicinity of the parents. In this particular optimization problem, each individual evaluation typically required ∼ 90 s what amounted for about 7.5 h of CPU time to perform a single run of the GA with 300 individuals [see Fig. 4]. We used a cluster of 12 cores within a Grid infrastructure, which reduced the computation time to ∼ 40 minutes. The Grid protocols supporting the GA execution make infrastructure resizable according to the needs of the problem: number of executions, dimensionality of the search space, etc. (see Ref. [19] for details on the Grid).

Dual-pulse solitonic source optimization
As mentioned above, the desired output spectral channels in this work intend to cover OCT applications in the NIR II window, where transparency of the biological tissues increases and scattering decreases [47]. Because spectral bell-shaped pulses avoid spurious structures in OCT images [48], the bright optical solitons are very good candidates for OCT applications. Figure 3 shows the spectral evolutions (bottom) and output spectrograms (top) corresponding to the best individuals obtained by the GA strategy and fitness function, Eq. 4, described in the previous section. All output spectra shown in Fig. 3 present the two reddest solitonic pulses, ejected from the soliton fission, accurately centered in the predefined channels (λ c 1,2 ) delimited by the dashed lines (see Table I for parameter values associated to results in Fig. 3). In Figs. 3(a) and 3(d), the target spectral channels where chosen from Ref. [29] in order to illustrate the solution for a λ (nm) (d-f) Spectral evolutions along the fiber associated to (a-c), respectively, retrieved from the optimal input pulse parameters corresponding to three different pairs of channels, λ c 1,2 . Input pulse parameters are given in Table I 3(c,f), demonstrate the tunability of such source, keeping λ c 1 − λ c 2 fixed to 100 nm without replacing the PCF but merely adjusting the input pulse parameters. We checked by benchmarks that several runs of the GA with fixed λ c 1,2 provided systematically very similar optimal results and therefore only one is shown here for each different case. Regarding OCT applications, another important aspect of the source presented here is that the fs-SC dynamics typically exhibits a very high coherence and negligible shot-to-shot fluctuations [3], known to be detrimental for OCT [35]. Moreover, the two output solitonic pulses (S 1,2 in Fig. 3) constituting the proposed OCT light source, provide a decent resolution l c ≡ 2 ln 2λ s /[π∆λ s,FW HM ] [49] of ∼ 10 µm for the two solitons, S 1,2 . Table 1. Parameters associated to the best individuals found by the GA, shown in Fig. 3.

Optimal pulse parameters
Spectral bands Resolution Fitness Shown in   curve in Fig. 4(b)]. However, the scattering ability of GOs often results in finding slightly better individuals in nearby regions presenting smaller agglomeration. An important reason for the convergence of our GA towards the optimal solutions is the fact that the operatorM is given a lower probability of action thanX (probabilities are 1/3 and 0.95 respectively, see previous section). This combination gives both a good diversity and probability to conserve the properties of the best individuals during the execution of the GA.

Conclusions
We presented an efficient optimization procedure based on GAs deployed in the Grid platform, providing faster results and potential scalability of the computational resources. The optimization provides the optimum input pulse parameters required to control the SC dynamics in a way that the first two ejected Raman solitons are centered at two pre-defined wavelengths. The results are shown to be of interest for practical OCT applications in the NIR II region where dual frequency, pulsed sources enable in vivo imaging, and avoid spurious results.