Ultrashort Pulse Simulations with the RP Fiber Power Software

The essential concept for simulating ultrashort pulse propagation in RP Fiber Power differs from that of our earlier product RP ProPulse. We do not first describe our optical system – the racetrack through which we later on want to send our pulses – to the software. Instead, we use certain functions of the script language for defining or loading an initial pulse and later on for sending that pulse through optical elements. Additional functions can be used for retrieving all sorts of properties of the “current pulse”, meaning the pulse as it is at the moment (after having seen certain elements).

The advantage of that approach is that it is most flexible. The mentioned functions can be used within any expressions, and we can make use of the powerful control structures available for such expressions: if-then statements, while-do loops, for loops, repeat-until and the like. You can use such expressions in many ways – for example, for various initializations, but also during the generation of any graphical diagram. Such flexibility turns out to be essential for any real research or development work.

In most cases, what is of main interest is what pulses you get in the steady state of the mode-locked laser. Usually, one cannot directly calculate that steady state, since it results from complicated interactions involving optical gain and losses (possibly wavelength-dependent), chromatic dispersion, fiber nonlinearities, saturable absorption, etc. Therefore, in most cases one just simulates the evolution of some initial pulse within some number of resonator round trips, hoping that the system will closely reach the steady state within a reasonable number of round trips. Well, there are cases where there is no such steady state; there are unstable lasers, the study of which can of course also be scientifically interesting. In some cases, one also likes to study the details of the start-up process, but let us now focus on the steady state.

The details of the initial pulse are usually not critical; you just don't want to be far away from the steady state in order to reach that quickly. You would just roughly guess parameters like pulse energy and duration, and generate a Gaussian unchirped start pulse, for example:

; Parameters:
l_s := 1060 nm  { center wavelength }
T_range := 25 ps  { width of temporal range }
N_t := 2^10  { number of grid points }
dt := T_range / N_t  { temporal resolution }
E0 := 10 nJ  { initial pulse energy }
tau0 := 1 ps  { initial pulse duration }
chirp0 := 0 GHz / ps  { initial chirp }
 
calc  { evaluate the following composite expression }
  begin
  set_pulse_grid(T_range, N_t, l_s);  { define the pulse grid }
  startpulse_G(E0, tau0, chirp0);  { generate the start pulse }
  end

The pulse grid, essentially defined by the width of the temporal range, the number of numerical points and the center wavelength, must be large enough to well accommodate the pulse at all stages. Also, the temporal resolution must be high enough to accommodate the full pulse bandwidth. For many bulk laser simulations, 256 points are fully sufficient, but for strongly chirped pulses, as often occur in fiber lasers, 1024 points or even more can be necessary. (For supercontinuum generation, one may require even tens of thousands.)

At some point, you have to define the involved fibers. Essentially, for each one you define a model as used for continuous-wave simulations, and in addition you assign some properties to it which are needed for ultrashort pulse propagation: essentially, the chromatic dispersion and the nonlinear index – in some cases, some more details concerning stimulated Raman scattering (not to be considered here). The following is the used code, not including the definition of various parameters, for a model with an active and one passive fiber in a linear resonator:

calc  { define the fibers }
  begin
    { active fiber: }
    set_device(1);
    set_fiber(L_active, 20, 'Yb');
    add_ring(r_co, N_Yb);
    pump := addinputchannel(P_pump_in, l_p, 'I_p', 0, forward);
    signal_active_fw := addinputchannel(0, l_s, 'I_s', 0, forward);
    signal_active_bw := addinputchannel(0, l_s, 'I_s', 0, backward);
    set_GVD(GVD_af);
    set_n2(n2_f);  { nonlinear index }
    finish_fiber();
    { passive fiber: }
    set_device(2);
    set_fiber(L_passive, 5, '-');
    signal_passive_fw := addinputchannel(0, l_s, 'I_s', 0, forward);
    signal_passive_bw := addinputchannel(0, l_s, 'I_s', 0, backward);
    set_GVD(GVD_pf);
    set_n2(n2_f);  { nonlinear index }
    finish_fiber();
  end;

The start pulse is usually assumed to be the intracavity pulse just before it hits the output coupler. So the output coupler, causing some losses, will be the first optical element seen in a resonator round trip. After that, the circulating pulse may see the active fiber, one or more passive fibers, a saturable absorber and/or an optical modulator, and possibly other optical elements. It is often convenient to define a function which performs a complete resonator round trip:

DoResonatorRoundTrip() := 
  { Simulate one resonator round-trip of the pulse. }
  begin
    global allow all;
    pp_multiply_expr_f('r_FBG%(f)');
      { fiber Bragg grating as output coupler }
    pp_fiber(1, signal_active_fw);  { active fiber }
    pp_fiber(2, signal_passive_fw);  { passive fiber }
    pp_sat_abs(dR_S, tau_S, E_sat_S);  { SESAM }
    pp_loss(loss_S);  { loss in SESAM }
    pp_fiber(2, signal_passive_bw);  { passive fiber }
    pp_fiber(1, signal_active_bw);  { active fiber }
    pp_center(1);  { center the pulse }
    calc_dyn(0, T_rt, T_rt);
      { dynamical simulation of gain recovery }
  end

Note that our example is based on a linear resonator, where after reaching the other end (the SESM, a saturable absorber) we go all the way back to the initial position. Of course, various parameters used in our function have to be defined elsewhere.

The last function call simulates the pumping of the fiber during one round-trip time, replenishing the stored energy.

Thereafter, it is easy to simulate some number of round trips and display some parameters of the intracavity pulse after those:

calc
  for j := 1 to 1000 do
  begin
    DoResonatorRoundTrip();
    store_pulse(j);
  end
show "Energy:    ", E_p():d3:"J"
show "Duration:  ", tau_p():d3:"s"
show "Bandwidth: ", dl_p():d3:"m"

Here, we have also stored all the pulses, so that we can later recall them for additional calculations and for generating diagrams, and also we can inspect them with the interactive pulse display window.

There is actually a problem with the simple approach shown here: usually it takes a very large number of round trips until the fiber's gain has settled at the final level. Essentially, there are relaxation oscillations which disappear only after a huge number of round trips. If you want to study those, it's just fine. If you want to find a steady state quickly, however, you can apply some additional trick: with a few more lines of code (not shown here), you can in each round trip have the fiber's gain recalculated according to the pulse energy – just as if the gain could respond to that instantaneously. That way, you can effectively suppress relaxation oscillations and get the final result much more quickly. Another possibility would be to artificially increase the gain saturation by a large factor (e.g. 100) while decreasing the effective pulse repetition rate by the same factor; that would effectively accelerate the relaxation oscillations very strongly, hugely saving on computation time. Of course, users of our software get my detailed advice on such tricks.

Another issue is that you may not be sure how many round trips are required for reaching the steady state. Therefore, you may want to define a function which simulates such round trips until the essential pulse properties (energy, duration, bandwidth) do not change substantially anymore:

FindSteadyState(tol, N_max) := 
  { Find the steady state,
    defined by the time where the relative changes
    of E_out and tau are less than tol for 3 successive steps.
    Return 1 if steady state is found within N_max steps,
    0 otherwise. }
  begin
    var n_sc, N, E1, E2, tau1, tau2, df1, df2;
    N := 0;  { number of round trips }
    n_sc := 0;
      { number of round trips with small change of pulse parameters }
    E1 := E_p();  tau1 := tau_p();  df1 := df_p();
    repeat
      inc(N);
      DoResonatorRoundTrip();
      E2 := E_p();  tau2 := tau_p();  df2 := df_p();
      if abs(E2 / E1 - 1) <= tol
         and abs(tau2 / tau1 - 1) <= tol
         and abs(df2 / df1 - 1) <= tol
      then inc(n_sc)
      else n_sc := 0;
      E1 := E2;  tau1 := tau2;  df1 := df2;
    until n_sc = 3 or N >= N_max;
    (n_sc = 3);
  end

You can then call that function at any time to get the steady state calculated efficiently. It will return 1 if it worked, or 0 in case that the steady state could not be found within the given maximum number of round trips. You could easily modify the details of that strategy as required in your case.

I give you an example for the application of that function. Let us assume that we want to generate a diagram where we want to plot the steady-state pulse energy and duration as functions of the pump power:

"Variation of Pump Power"
 
x: 400, 600
"pump power (mW)", @x
y: 0, 20
frame
 
! for x := CS_x1 to CS_x2 step 20 do
  begin
    set_device(1);  { active fiber }
    set_P_in(pump, x * mW);
    if FindSteadyState(0.01, 1000) then
    begin
      point(x + i * E_p() / nJ, "R");
        { filled rectangle for pulse energy in nJ }
      point(x + i * tau_p() / ps, "t");
        { open triangle for pulse duration in ps }
    end;
  end

The example above has shown you how one can implement a relatively sophisticated simulation with not too much code. Here are some more example for possible scenarios:

• Consider an amplifier system fed by a continuous pulse train from a seed laser. Again, you can simulate how the pulse energy evolves towards its steady state. You could also study e.g. how it fluctuates if you have random time intervals between the pulses, or just some pulses missing, or a sequence of pulse bunches, or fluctuations of the pump power or of properties of the seed pulses, etc.
• A regenerative amplifier works such that one sometimes injects a seed pulse, then lets it evolve over some number of round trips, and then ejects it. After some time of pumping the amplifier crystal, one can perform the next amplification cycle. Here, one can easily implement a function which simulates a complete amplification cycle with the given number of resonator round trips. Again, one could use a function which simulates such amplification cycles until the steady state is reached. Some systems exhibit bifurcations and even chaos – see some online example case.

In the example above, you have seen how to send a pulse through different optical fibers, and also how to apply optical components with some optical losses or some saturable absorption. There are lots of other functions for applying other optical components:

• pp_dispersion() for applying chromatic dispersion up to 4th order
• pp_multiply_expr_f() for multiplying the frequency-domain amplitudes of the pulse with a frequency-dependent factor – for example, for applying arbitrary dispersion (phase shifts) or a bandpass filter; similar things can be done in the time domain, e.g. for optical modulators
• pp_prism_pair() and pp_grating_pair() for applying the chromatic dispersion from a prism pair or grating pair, respectively, with given parameters (such as spacing, lines per millimeter, etc.)
• pp_compress() for an automatically optimized dispersive compressor (up to 4th order)
• pp_noise() for adding some random noise
• pp_SPM() for an optical element with self-phase modulation
• pp_add_pulse() for coherently adding a pulse, e.g. in simulations for interferometers or of additive pulse mode locking

These components allow you to simulate all common kinds of setups for the generation or manipulation of ultrashort pulses, only (so far) not with parametric nonlinear gain.

All that has been developed with great care. To give you an example, I realized that functions like pp_multiply_expr_f() and pp_prism_pair() are often applied many times with the same parameters. I therefore implemented those such that the frequency-dependent factors, as obtained to the frequency-domain pulse amplitudes, are stored and can be reused in a later function call. That works even if multiple parameter sets are used in a simulation. Due to such tricks, which the user does not even have to think about, the performance of such computations can be very high.

I hope you have seen that it can be great fun to have such a powerful tool. A researcher can use it to relatively quickly develop rather sophisticated studies, and in industry you can rapidly analyze the behavior of systems in order to optimize them or understand certain peculiarities. The alternative would be working in blind flight mode – not very efficient in practice …

Of course, ultrashort pulse simulations are not a trivial thing, involving a lot of physics and also a lot of important practical aspects. Therefore, it is essential not only to have some good piece of software, but also to obtain powerful technical support. I deliver that personally, making sure that people can make rapid progress, and that way I regularly have very happy customers.