# Supporting Evolution of Safety Cases within
the Software Development Process

Ankit Agrawal  Department of Computer Science
Saint Louis University
Saint Louis, MO, USA
<EMAIL_ADDRESS>Jane Cleland-Huang  Department of Computer Science
University of Notre Dame
SouthBend, IN, USA
<EMAIL_ADDRESS>
# Leveraging Traceability to Integrate Safety Analysis Artifacts into the
Software Development Process

Ankit Agrawal  Department of Computer Science
Saint Louis University
Saint Louis, MO, USA
<EMAIL_ADDRESS>Jane Cleland-Huang  Department of Computer Science
University of Notre Dame
SouthBend, IN, USA
<EMAIL_ADDRESS>
###### Abstract

Safety-critical system’s failure or malfunction can cause loss of human lives
or damage to the physical environment; therefore, continuous safety assessment
is crucial for such systems. In many domains this includes the use of Safety
assurance cases (SACs) as a structured argument that the system is safe for
use. SACs can be challenging to maintain during system evolution due to the
disconnect between the safety analysis and system development process.
Further, safety analysts often lack domain knowledge and tool support to
evaluate the SAC. We propose a solution that leverages software traceability
to connect relevant system artifacts to safety analysis models, and then uses
these connections to visualize the change. We elicit design rationales for
system changes to help safety stakeholders analyze the impact of system
changes on safety. We present new traceability techniques for closer
integration of the safety analysis and system development process, and
illustrate the viability of our approach using examples from a cyber-physical
system that deploys Unmanned Aerial Vehicles for emergency response.

###### Index Terms:

Safety Case, Safety Analysis, Traceability

##  I Introduction

Safety-critical systems are systems whose failure could result in loss of
life, significant damage to the environment, or significant financial loss
[11]. Such systems must be developed systematically and rigorously. Given a
set of requirements describing the system’s functionality, we need to assure
that associated hazards have been identified and appropriately addressed,
typically using techniques such as Fault-Tree Analysis (FTA) and Failure-Mode
Effect and Criticality Analysis (FMECA). Beyond these techniques, it is
increasingly common for organizations to construct claim-based safety
arguments [12] in the form of a Safety Assurance Case (SAC). A SAC decomposes
high-level safety goals or claims into layers of arguments supported by safety
evidence such as test-cases logs, simulation results, or formal proofs [3],
often using either the Claims-Arguments-Evidence notation [16] or the Goal
Structuring Notation [13].

SACs are recommended, or even required, in many safety critical domains
(e.g.,[19]); however, in a study by Cheng et al., safety experts reported that
there are ‘no effective mechanisms for managing change’ for a SAC, and that
‘SAC creation and maintenance has not been fully integrated into the software
development process’ [5]. System change analysis becomes difficult because of
insufficient domain-knowledge among safety stakeholders and the sheer
complexity involved in safety-critical products. During change analysis, the
safety stakeholders seek answers to questions such as (1) why the system has
changed, (2) what risk does this change mitigates, and (3) how this change can
impact safety[10]. Further, safety stakeholders typically design and maintain
SACs while the development team produces artifacts, such as mitigating
requirements and test results, on which the SAC’s safety arguments depend.
These distinct roles also create a gap between system development and SAC
maintenance processes. Therefore, it is crucial not only to establish
traceable links between various safety artifacts such as Fault Tree, SAC, and
development artifacts such as Design Decisions and code, but to keep safety
analysts informed of how changes in development artifacts may affect safety by
documenting the rationales for changes to these development artifacts.

In this paper, we propose a solution to establish traceability links between
various safety artifacts such as SACs and Fault Trees, and supplementing the
SACs with rationales for changes in the development artifacts. Our goal is to
improve maintainability of SACs and keep safety analysts informed of the
rationale for changes in the development artifacts as the software evolves.
Our solution first utilizes Safety Artifact Forest Analysis (SAFA) [2, 17, 8],
which detects changes in the software development artifacts between the two
versions of the underlying system and automatically generates vizualizations
enabling analyst to easily navigate through changes. We establish traceable
links between these auto-generated vizualizations and safety artifacts such as
FTAs and SACs. Secondly, we discuss strategies to capture the reasons behind
changes in software development artifacts and maintain change rationale
information as part of the development artifacts. Finally, we demonstrate how
this comprehensive traceability across multiple safety artifacts, supported by
rationales for changes in the development artifacts, enables us to analyze the
impact of changes on safety as the system evolves. To illustrate our approach,
we provide examples from the DroneResponse system [7, 1, 6], which utilizes
Unmanned Aerial Vehicles (UAVs) to support emergency response. Finally, we
outline the open challenges ahead in a preliminary roadmap.

##  II System Artifacts Traceability

Figure 1: Software artifacts, fault tree nodes, and evidence nodes from the
SAC are connected via trace links to integrate system development and safety
processes and to provide support for managing the impact of change upon system
safety.

The SAFA framework retrieves artifacts from project repositories such as
DOORS, Jira, and Github. Given a root node, such as a system level
requirement, it constructs a vertical slice through the system according to
the trace links defined in a Traceability Information Model (TIM). SAFA refers
to such a tree as an Artifact Tree (AT). SAFA can compare a current version of
the AT against an earlier baseline version to produce a Delta Tree (DT) which
visualizes changes in the system. The right lane of Figure 1 shows a partial
DT. Additional examples are provided in our prior work [2].

SAFA detects additions (green), deletions (red), and modifications (blue) for
requirements, design, code, tests, operating context, environmental
assumptions, and other system artifacts. The delta tree visualization helps
project stakeholders to identify changes and to investigate their impact on
system safety. It highlights these changes and recommends areas in which
inspection is needed. Safety experts expressed the importance of providing
rationales for each change [2]. While SAFA currently uses basic static
analysis techniques to explain refactoring changes; richer rationales are
needed that describe changes in requirements and design, as well as
modifications to the code.

We illustrate the proposed solution with examples from publicly available
requirements dataset for DroneResponsecyber-physical system that deploys
cohorts of UAVs to support emergency response mission such as search-and-
rescue [9]. The use of on-board intelligence by UAVs for autonomous navigation
through airspace is one of the key requirements of multi-UAV autonomous
systems. However, the US Federal Aviation Authority (FAA) regulates airspace
usage and defines special-use or restricted-airspace where UAVs are not
allowed to fly. Autonomous UAVs entering the restricted-airspace could cause
accidents with commercial flights, military operations, or medivac deliveries.
Therefore, we consider flights into prohibited space to be a severe
operational risk and have carefully designed mitigations into our system.

The safety requirement UAV-1387 in Figure 1 states that “When the UAV is in
autonomous mode it shall fetch restricted airspace data from the LAANC
system.” In the previous version of the system, this was partially addressed
by requiring the UAV to continuously check for airspace information while in
autonomous mode (Design Definition - UAV 1388). However, in the new system,
this requirement was replaced by conducting a more economical check when new
flight paths are planned (Design Definition - UAV 1413). The delta tree,
produced by SAFA, clearly highlights this design change, showing the
replacement of Design Definition UAV-1388 and its associated code in red, and
the inclusion of Design Definition UAV-1413 and its associated code in green.

##  III Integrating System Artifacts with Safety Assets

When developing a safety-critical system, a preliminary hazard analysis (PHA)
is performed [14] to identify high-level hazards that represent undesirable
states of the system. Each of these hazards is then explored through an
associated Fault Tree (FT) [18] or a FMECA model [15]. In this paper, we
illustrate our approach using FTs; however, our techniques are also applicable
to FMECAs. Fault Tree Analysis (FTA) is a top-down approach that starts by
analyzing the high-level risk and then uses boolean logic to depict a chain of
events causing system-level risk. The middle lane of Figure 1 shows a partial
FT for risks associated with UAV flights in restricted airspace.

FTs and FMECAs are typically used as part of a SAC’s argumentation structure
to show that a specific fault has been sufficiently mitigated. Therefore a
link should be established from a FT to its relevant argument in the SAC. An
example of a SAC argument using GSN is depicted in the left lane of Figure 1.
This SAC uses a complete fault tree as evidence to support the claim that
hazards due to FAA restrictions on airspace are mitigated. To maintain
horizontal traceability (depicted by horizontal dashed arrows in Figure 1), an
explicit link is established between an evidence node in the SAC and the root
node of the FT. In turn, links are established from multiple nodes in the FT
to sub-trees of system artifacts. For example, the intermediate fault node of
_‘drone autonomously navigates into the restricted airspace’_ is linked to an
acceptance test as verification that the fault has been mitigated. Similarly,
one of the contributing basic faults _‘Drone operates on stale restricted-
airspace data’_ is linked to the previously discussed safety requirement
(UAV-1387).

In our solution, we establish traceability links between system artifacts,
safety assets (e.g., FTs and FMECAs), and SACs to propagate changes back-and-
forth between safety assets and system artifacts. Figure 1 illustrates that a
single change in the design of the system could trigger and propagate
notifications and warnings across the linked safety artifacts (yellow nodes).

##  IV Capturing Rationales

To support Safety Analysts in analyzing the impact of changes in system
artifacts on the overall safety of the system, our approach strategically
captures rationales from developers and other project stakeholders as they
make changes that impact artifacts linked directly or indirectly to an FT,
FMECA, or other safety assets. For example, in the FT depicted in the middle
lane of Figure 1, we observe that the leaf node describing the basic fault
_‘Drone operates on stale restricted airspace data’_ is mitigated through
safety requirement UAV-1387, which states that _‘When the UAV is in autonomous
mode it shall fetch restricted airspace data from the LAANC system.’_ As
indicated by the red nodes, the original design satisfied this requirement
through continuously fetching restricted airspace data during flight, while
the current version (shown in green) replaces this functionality with a single
fetch each time a new flight path is planned. A safety analyst will need to
determine whether this change adversely impacts safety.

Burge et al suggested capturing reasons, alternative options, and arguments to
provide rationales for design decisions [4]. We therefore elicit a rationale
for any design decision that links directly to an FTA. In the case of
DroneResponse, all requirements and design decisions are captured in Jira, and
Figure 2 shows possible instrumentation of the Jira environment to elicit
design rationales when the original design requirement (UAV-1388) is replaced
by a new one (UAV-1413).

Figure 2: User interface for capturing design decision details, including
reasons, alternatives, and arguments. Figure 3: Eclipse Plugin User Interface
for capturing code-level change rationale

At the code level we capture change rationales. These are akin to commit
messages, but at the granularity of each modified class instead of an entire
change set. In our example, the MonitorAirspace.java file is replaced by a new
optimized service OnDemandAirspace.java that fetches data on-demand. Figure 3
shows the user interface of a prototype IDE plugin that could be used to
capture granular details of the change. The plugin not only elicits a
justification and explanation of the change from the developer, but also
visualizes the contribution that the class makes to an FTA (e.g., Fig. 1) and
provides links to prior commit messages and previous change justifications.

The rationales captured as part of the change process provide crucial domain
knowledge to conduct a thorough analysis of system changes on safety.
Therefore, the domain knowledge acquired from these rationales, alongside
visual representations of change impact, as depiceted in Figure 1, can aid
safety analysts to determine (1) whether current changes impact safety or not,
(2) whether additional mitigations are needed, and/or (3) whether the FT,
FMECA, and/or the SAC need to be updated to reflect new hazards or new safety
arguments. Decisions made by safety analysts can then be propagated back to
the development team in order to close the loop.

##  V Future Challenges

  * •

Knowledge Management: As previously reported, the innate complexity of safety-
critical systems means that experts from diverse disciplines are needed to
construct and maintain a SAC [5]. While we integrate basic rationale capture
into SAFA, open questions include (1) What types of domain knowledge are
needed by safety analysts to evaluate and/or construct a SAC? (2) How and when
should this information be collected? (3) How can it be effectively used to
support Safety Analysts?

  * •

Intelligent analysis of change:  Changes are introduced into a system at many
different levels to accommodate changes in the environment, introduce
functional enhancements, improve system qualities such as performance,
reliability, or maintainability, or to correct errors. From a safety
perspective we need to differentiate between harmful and non-harmful errors.
Our current approach visualizes all change, captures stakeholders’ change
rationales, and provides basic explanations for change. However, future
systems should be able to leverage AI solutions to (1) analyze individual and
composite changes in order to identify and explain patterns of change, (2)
differentiate between harmless changes and those with potential safety impact,
and (3) recommend remedial actions when safety is impacted.

  * •

Tool Supported Integrated Environments Our proposed approach requires trace
links to be established across heterogeneous artifacts stored in diverse tools
and repositories. Safety experts have previously reported that the lack of
tool support and clear guidance make SAC creation and maintenance challenging
[5]. Open challenges therefore include (1) defining best practices for
creating end-to-end traceability across safety assets and development
artifacts, (2) providing integrated tool-supported environments for retrieving
diverse artifacts and establishing effective traceability that supports safety
analysis, and (3) developing interactive visualization tools that display
SACs, FTAs, FMECAs, rationales, and artifacts in ways that provide appropriate
support for different types of users such as safety analysts and developers.

This paper has presented a preview of our current research in addressing the
disconnect between safety analysis and the software development process. As
proposed in this paper, we are developing plugins to capture change rationales
at the design and code level, and extending SAFA to show rationales in Delta
view and link FTAs, FMECAs, and SACs in it.

## References

  * [1] A. Agrawal, S. J. Abraham, B. Burger, C. Christine, L. Fraser, J. M. Hoeksema, S. Hwang, E. Travnik, S. Kumar, W. Scheirer, et al.  The next generation of human-drone partnerships: Co-designing an emergency response system.  In Proceedings of the 2020 CHI Conference on Human Factors in Computing Systems, pages 1–13, 2020.
  * [2] A. Agrawal, S. Khoshmanesh, M. Vierhauser, M. Rahimi, J. Cleland-Huang, and R. Lutz.  Leveraging artifact trees to evolve and reuse safety cases.  In 2019 IEEE/ACM 41st International Conference on Software Engineering (ICSE), pages 1222–1233. IEEE, 2019.
  * [3] P. Bishop and R. Bloomfield.  A Methodology for Safety Case Development.  In Ind. Perspect. Safety-critical Syst., volume 19, pages 194–203. Springer London, 1998.
  * [4] J. E. Burge and D. C. Brown.  Software engineering using rationale.  Journal of Systems and Software, 81(3):395–413, 2008.
  * [5] J. Cheng, M. Goodrum, R. Metoyer, and J. Cleland-Huang.  How do practitioners perceive assurance cases in safety-critical software systems?  In Cooperative and Human Aspects of Software Engineering, pages 57–60. ACM, 2018.
  * [6] J. Cleland-Huang and A. Agrawal.  Human-drone interactions with semi-autonomous cohorts of collaborating drones.  arXiv preprint arXiv:2010.04101, 2020.
  * [7] J. Cleland-Huang, A. Agrawal, M. N. A. Islam, E. Tsai, M. Van Speybroeck, and M. Vierhauser.  Requirements-driven configuration of emergency response missions with small aerial vehicles.  In Proceedings of the 24th ACM Conference on Systems and Software Product Line: Volume A-Volume A, pages 1–12, 2020.
  * [8] J. Cleland-Huang, A. Agrawal, M. Vierhauser, and C. Mayr-Dorn.  Visualizing change in agile safety-critical systems.  IEEE Softw., 38(3):43–51, 2021.
  * [9] J. Cleland-Huang, M. Vierhauser, and S. Bayley.  Dronology: an incubator for cyber-physical systems research.  In A. Zisman and S. Apel, editors, Proceedings of the 40th International Conference on Software Engineering: New Ideas and Emerging Results, ICSE (NIER) 2018, Gothenburg, Sweden, May 27 - June 03, 2018, pages 109–112. ACM, 2018.
  * [10] M. Gleirscher and C. Carlan.  Arguing from hazard analysis in safety cases: a modular argument pattern.  In 2017 IEEE 18th International Symposium on High Assurance Systems Engineering (HASE), pages 53–60. IEEE, 2017.
  * [11] W. S. Greenwell, E. A. Strunk, and J. C. Knight.  Failure analysis and the safety-case lifecycle.  In Human Error, Safety and Systems Development: IFIP 18th World Computer Congress TC13/WC13. 5 7th Working Conference on Human Error, Safety and Systems Development 22–27 August 2004 Toulouse, France, pages 163–176. Springer, 2004.
  * [12] R. Hawkins, I. Habli, T. Kelly, and J. McDermid.  Assurance cases and prescriptive software safety certification: A comparative study.  Saf. Sci., 59:55–71, 2013.
  * [13] T. Kelly and R. Weaver.  The Goal Structuring Notation – A Safety Argument Notation.  In Proc. Dependable Syst. Networks 2004 Work. Assur. Cases, 2004\.
  * [14] N. G. Leveson.  Safeware, System Safety and Computers.  Addison Wesley, 1995.
  * [15] R. R. Lutz and R. M. Woodhouse.  Requirements analysis using forward and backward search.  Ann. Software Eng., 3:459–475, 1997.
  * [16] K. Netkachova, O. Netkachov, and R. Bloomfield.  Tool support for assurance case building blocks: Providing a helping hand with cae.  In Computer Safety, Reliability, and Security: SAFECOMP 2015 Workshops, ASSURE, DECSoS. ISSE, ReSA4CI, and SASSUR, Delft, The Netherlands, September 22, 2015, Proceedings 34, pages 62–71. Springer, 2015.
  * [17] A. D. Rodriguez, T. Newman, K. R. Dearstyne, and J. Cleland-Huang.  SAFA: A tool for supporting safety analysis in evolving software systems.  In 37th IEEE/ACM International Conference on Automated Software Engineering, ASE 2022, Rochester, MI, USA, October 10-14, 2022, pages 165:1–165:4. ACM, 2022.
  * [18] N. R. Storey.  Safety Critical Computer Systems.  Addison-Wesley Longman Publishing Co., Inc., Boston, MA, USA, 1996.
  * [19] U.S. Food and Drug Administration.  Infusion Pumps Total Product Life Cycle: Guidance for Industry and FDA Staff, 2014.


# Massive star feedback in the Magellanic Clouds and the tidal Bridge

Varsha Ramachandran 1 1Zentrum für Astronomie der Universität Heidelberg,
Astronomisches Rechen-Institut, Mönchhofstr. 12-14, 69120 Heidelberg

(2022)

###### Abstract

Massive stars have far-reaching feedback effects that alter the surrounding
environment on local, global, and cosmic scales. Spectral analyses of massive
stars with adequate stellar-atmosphere models are important to study massive
star feedback in detail. We discuss the most recent UV and optical studies of
massive metal-poor stars, including those with metallicities ranging from half
to one twentieth of solar, connected with large-scale ISM structures in the
Magellanic Clouds and the tidal Magellanic Bridge. We present ionizing fluxes
from massive stars with low metallicity along with mechanical energy, and we
further compare these to the observed energetics in the ISM. The results give
hints on the leakage of hot gas and ionizing photons in the Magellanic Clouds.
The paper outlines feedback from individual massive stars to population-level
collective feedback, the significance of various feedback mechanisms
(radiation, wind, supernova), and the influence by the physical conditions of
the ISM.

###### keywords:

Massive stars, stellar feedback, low metallicity, etc.

††journal: Astronomy in Focus, Focus Meeting 4††editors: José Espinosa, ed

##  1 Introduction

Stellar feedback is one of the largest uncertainties in star and galaxy
formation. Massive stars are of great interest because of their far-reaching
feedback effects that alter the surrounding environment on local, global, and
cosmic scales. The combined feedback in massive young stellar clusters leads
to the formation of superbubbles that can drive galactic winds and outflows,
which have been frequently observed in local as well as distant galaxies. The
Magellanic Clouds and the Bridge offer an outstanding opportunity to
investigate low metallicity massive stars and study feedback under conditions
typical for the vast majority of dwarf galaxies.

##  2 Massive star feedback in the Magellanic Clouds

We carried out spectroscopic observations of $\sim\\!500$ OB stars in the
Magellanic Clouds using VLT-FLAMES. When available, the optical spectroscopy
was complemented by UV spectra from the HST, IUE, and FUSE archives. The two
representative young stellar populations that have been studied are associated
with the superbubble N 206 in the Large Magellanic Cloud (LMC) and with the
supergiant shell SMC-SGS 1 in the Wing of the Small Magellanic Cloud (SMC),
respectively. We performed spectroscopic analyses of the massive stars using
the non-LTE Potsdam Wolf-Rayet (PoWR) model atmosphere code. We estimated the
stellar, wind, and feedback parameters of the individual massive stars.

The total energy feedback from OB stars, WR stars and supernovae in the N 206
and SGS complexes are compared in Table 1. Since the winds of OB stars in the
LMC are much stronger than in the SMC, wind feedback in N 206 is dominated by
massive Of stars (Ramachandran et al., 2018a, b). The stellar wind feedback
from the two WR stars in N 206 is comparable to that from all young OB stars
together. The situation is completely different in the SMC supergiant shell,
where the ionizing and mechanical luminosity is dominated by one WO star
(Ramachandran et al., 2019). The observed H$\alpha$ emission reflects almost
half of the ionizing flux provided by the massive stars – the other half of
the LyC photons obviously escapes. The X-ray superbubble in the LMC is equally
powered by stellar winds and supernovae. In contrast, the overall energy
feedback of the supergiant shell in the SMC is dominated by supernovae,
implying that hundreds of OB stars do not have much impact on the feedback
before their final explosion. Thus, at low metallicities, feedback is chiefly
governed by supernova explosions and only few very massive stars. The
comparison of the total stellar input with X-ray, radio and H$\alpha$
observations shows that only a fraction the input energy accumulated over time
is currently still present in these regions. The rest might have escaped or
leaked out of the complex. Our result has significant importance for feedback
in low-metallicity galaxies, where we can now neglect the contribution from OB
stars. However, at LMC like metallicities, young OB stars have a significant
contribution in terms of ionizing and mechanical feedback. In that case,
neglecting their contribution and considering only supernovae is not
justified. Concluding, the metallicity decides whether the stellar winds or
supernovae become the key agents of feedback.

Total $E_{\mathrm{mec}}$ | N206-LMC | SMC-SGS 1
---|---|---
OB |  $10^{52}$ (43%) |  $6\times 10^{49}$ (0.5%)
WR |  $10^{51}$ (9%) |  $2\times 10^{51}$ (16%)
SNe |  $1.2\times 10^{52}$ (48%) |  $10^{52}$ (83%)
Observed |  $4.7\times 10^{51}$ (20%) |  $6\times 10^{51}$ (50%)
Table 1: Total mechanical energy feedback

Figure 1: Chemical abundance in Bridge O stars

##  3 Metallicity and ionization of the Bridge

The Magellanic Bridge, stretching between SMC and LMC, is the nearest tidally
stripped intergalactic environment. The Bridge has a significantly low average
metallicity than SMC. For the first time we discovered three massive O stars
in the Bridge using VL/FLAMES spectra (Ramachandran et al., 2021). We analyze
the spectra of each star using the PoWR models, providing stellar parameters,
ionizing photon fluxes, and surface abundances. The ages of the newly
discovered O stars suggest that star formation in the Bridge is ongoing. The
multi-epoch spectra indicate that all three O stars are binaries. Despite
their spatial proximity to one another, these O stars are chemically distinct.
MBO1 is a fast-rotating giant with nearly LMC-like abundances. The other two
are main-sequence stars that rotate extremely slowly and are strongly metal
depleted (Fig 1). Among this MBO2 the most nitrogen-poor O star known to date.
Taking into account the previous analyses of B stars in the Bridge, we
interpret the various metal abundances as the signature of a chemically
inhomogeneous ISM, suggesting that the Bridge gas might have accreted during
multiple episodes of tidal interaction between the Clouds. Attributing the
lowest derived metal content to the primordial gas, the time of the initial
formation of the Bridge may date back several billion years. Using the Gaia
and Galex color-magnitude diagrams, we roughly estimate the total number of O
stars in the Bridge and their total ionizing radiation. Comparing this with
the energetics of the diffuse ISM, we find that the contribution of the hot
stars to the ionizing radiation field in the Bridge is less than 10% and
conclude that the main sources of ionizing photons are leaks from the LMC and
SMC. This provide a lower limit for the fraction of ionizing radiation that
escapes from these two dwarf galaxies.

## References

  * Ramachandran et al. (2018a) Ramachandran, V., Hainich, R., Hamann, W.-R., et al. 2018a, A&A, 609, A7
  * Ramachandran et al. (2018b) Ramachandran, V., Hamann, W. R., Hainich, R., et al. 2018b, A&A, 615, A40
  * Ramachandran et al. (2019) Ramachandran, V., Hamann, W. R., Oskinova, L. M., et al. 2019, A&A, 625, A104
  * Ramachandran et al. (2021) Ramachandran, V., Oskinova, L. M., & Hamann, W. R. 2021, A&A, 646, A16


# Internal nonlinear transmission in an Yb mode-locked fiber laser through
bifurcations

Cecília L. A. V. Campos, Lucas B. A. Mélo,
Hugo L. D. de S. Cavalcante, Lúcio H. Acioli, and
Marcio H. G. de Miranda
Departamento de Física, Universidade Federal de Pernambuco
50740-560 Recife, PE - Brazil
Departamento de Informática, Centro de Informática,
Universidade Federal da Paraíba
58051-900 João Pessoa, PB - Brazil
<EMAIL_ADDRESS>
Published in Optics Communications, 26 10, 13686 (2019)

###### Abstract

Mode-locked fiber lasers are rich dynamical systems that may present several
different types of pulsed operating modes, depending on external control
parameters such as pump power. A systematic experimental characterization of
such regimes is a challenging problem. Here single pulse regimes of ytterbium
mode-locked fiber laser are explored and related to the nonlinear amplitude
modulation affecting the pulses inside the laser cavity. To understand the
full dynamics of this system, average power measurements (slow time scale)
have been performed and compared to simultaneous time-resolved measurements
(fast time scale). The average nonlinear transmission resulting from the
amplitude modulation is measured, allowing to relate the transitions between
dynamical regimes to the system’s nonlinearity.

##  1 Introduction

The interest in ytterbium-doped mode-locked fiber lasers (MLFL) has increased
over the last decades because, beyond been a relatively cheap system, it has
important technological applications in the fields of physics, engineering,
biology and medicine. It is also very interesting from the dynamical systems
point of view because it presents a dense forest of different phenomena to be
explored such as: deterministic chaos [1], rogue waves [2, 3], soliton
explosion [4], dissipative soliton molecules [5], harmonic mode-locking [6]
and multi-pulsing regime [7, 8, 9].

The Yb MLFL is quite flexible in the sense that the system operates in the
pulsed regime under very different conditions. One is allowed to externally
control two basic properties which determine the operating regime: (i) group
velocity dispersion (GVD) $\beta_{2}$, and (ii) the cavity nonlinear amplitude
modulation (nonlinear loss or saturable absorber). The parameter
$\left.\beta_{2}\right|_{\text{fiber}}$ of most optical fibers is positive at
$\lambda\approx 1~{}\mu$m, which is the Yb lasers emission. For Yb laser
emission most optical fibers present parameter
$\left.\beta_{2}\right|_{\text{fiber}}>0$. Even so, using a grating pair, one
can set the overall dispersion in the cavity at will. To obtain a stable mode-
locked pulse train, the GVD is usually set to $\beta_{2}\lesssim 0$, but short
pulse generation from Yb fiber lasers with overall normal dispersion,
$\beta_{2}>0$ regime has been demonstrated in references [10, 11].

The presence of nonlinear amplitude modulation inside the laser cavity is a
fundamental requirement to obtain mode-locked operation [12, 13], and in many
MLFL the most common choice to obtain intensity-dependent losses is based on
the nonlinear polarization rotation (NPR) [14], due to the nonresonant
$\chi^{(3)}$ of the optical fiber. This emulates a fast saturable absorber,
which has been extensively studied in the operation of pulsed femtosecond
lasers. In fiber lasers this is achieved by a combination of waveplates and
polarizers to obtain a nonlinear loss mechanism of the form

$L(P)=L_{0}-\gamma\,P\,,$ (1)

where $P$ is the instantaneous power, $L_{0}$ is the overall linear cavity
loss and $\gamma$ is the nonlinear loss coefficient. Suitable settings of
waveplates allow one to control both $L_{0}$ and $\gamma$. The inclusion of
higher order terms in Eq. (1) is essential for obtaining stable mode-locked
operation [16], and also known to be important to achieve different dynamical
regimes [15].

For a laser operating in the mode-locked regime due to NPR, the losses inside
the cavity are directly related to the nonlinear transmission occurring in a
polarizing element, such as a polarizing beamsplitter. Notice that the element
itself behaves linearly, with the nonlinearity originated in the fiber before
it.

The nonlinearities present in MLFL provide a fertile ground for experimental
exploration, and in order to understand their origin, one option has been to
rely on theoretical models based on the complex Ginzburg-Landau equation
(CGLE) or the nonlinear Schrödinger equation. This approach has been used to
study pulsating and chaotic solitons, period doubling, and pulse coexistence
in mode-locked lasers [17]. Wei et al. [18] have exploited the nonlinear
amplitude modulation to understand the physics behind the stochasticity,
hysteresis and multistable regimes. Komarov et al. [8] propose to study a
multistable Yb mode-locked laser focusing on the roles of group velocity
dispersion and the gain saturation.

A different approach was proposed by Li et al. [7], who have investigated the
importance of the nonlinear transmission to explain the transition between
single pulse and multi-pulse regimes within the fundamental repetition period
of a mode-locked laser. Our interest here is to experimentally characterize
different dynamical regimes, beyond stable mode-locked operation of an Yb
MLFL, and their relation to the nonlinear cavity transmission. This is done by
measuring average power values (slow time scale) over the whole pulse train
and comparing to time-resolved measurements (fast time scale). We have
discovered that important details about the nonlinear dynamics happening at
the fast time scale can be revealed by looking at the slow time scale
parameters. This is interesting and worthy of attention because slow detectors
are much more accurate and accessible than fast ones, so we suggest that the
nonlinear transmission curve can be an alternative for monitoring stability
and dynamical transitions of systems similar to the one presented here.

In this work, we have measured the nonlinear transmission curve for pulse
trains in an Yb MLFL, presenting a clear experimental signature of dynamical
regime transitions. To measure the nonlinear transmission curve, we have
monitored laser beams inside and outside of the laser cavity, which we will
refer as intracavity and extracavity, respectively. We have also measured the
average laser power, RF and optical spectra, and the temporal series. The
average laser power intracavity and extracavity presented a counter intuitive
tendency, where it does not grow linearly with the pump current but has a
completely sharp variation in slope due to the transition between dynamical
regimes. When both powers are added, they show a complementary behavior,
giving a linear growth.

##  2 Experimental setup

The experiment was performed using a homemade ytterbium-doped MLFL, similar to
the one used in Ref. [1]. A schematic diagram of the experimental setup and
data acquisition system can be seen in Fig. 1(a). The pump source is a diode
laser operating at 976 nm, with maximum power of 720 mW. The gain medium is a
highly Yb-doped single-mode fiber CorActive Yb214, with absorption coefficient
of 1348 dB/m at 976 nm and length of 22 cm. All the other fibers are standard
single-mode HI 1060. The pump power is coupled to the gain medium through a
wavelength division multiplexing (WDM), presenting a coupling efficiency of
89$\%$. A pair of GRIN collimators couple the beam to and from the free-space
section and provides optical feedback. The coupling efficiency back to the WDM
is 35$\%$. The intracavity dispersion is controlled by a grating pair with 600
grooves/mm, each grating presenting a efficiency of 65$\%$ per pass. An
optical isolator guarantees the unidirectionality of the signal. The laser is
passively mode-locked via NPR technique [14]. This technique is based on the
intensity-dependent evolution that a polarization state can exhibit during its
propagation in the fibers, provided that its intensity is sufficiently high.
The introduction of a polarizer into the cavity will cause an intensity
dependent loss, emulating a fast saturable absorber and allowing mode-locking
regime to occur. To adjust the polarization and achieve mode-locked operation,
a set of waveplates and a polarizing beamsplitter (PBS) are included (Fig.
1(a)). It is important to note that the PBS is also the output coupler of the
laser, meaning that the intracavity and the extracavity beams are linearly
polarized and orthogonal. The waveplates and PBS losses are very small and
will be neglected in the following discussion. When mode-locked, our laser
presents a train of pulses with repetition rate of 125 MHz, optical spectrum
centered in 1025 nm and temporal width of approximately 200 fs [1].

Figure 1: (a) Schematic representation of the experimental setup: diode laser
(DL), wavelength division multiplexing (WDM), ytterbium doped fiber (YDF),
quarter-wave plate ($\lambda$/4), half-wave plate ($\lambda$/2), polarizing
beamsplitter (PBS), grating pair (GP), Faraday isolator (FI), GRIN collimators
(C), power meter’s detectors (D1 and D4), 1 GHz photodetectors (D2 and D5) and
pick-up fiber for the OSA (D3). (b) Simple diagram of the laser beams in the
PBS.

The experiment consists in scanning the diode pump laser current
($I_{\text{pump}}$) from 525 mA to 969 mA, in steps of 4 mA, keeping the
waveplates settings and dispersion fixed. The calibration from
$I_{\text{pump}}$ to optical power is 0.68 mW/mA. Data acquisition was
automated and controlled by a custom LabVIEW program, involving two
photodetectors (1 GHz bandwidth), two RF spectrum analyzers (3 GHz bandwidth),
two power meters, one optical spectrum analyzer (OSA) and one oscilloscope (1
GHz bandwidth).

The time-resolved measurements are temporal series of pulses acquired by fast
photodetectors (1 GHz bandwidth). We denote this as a fast time scale
measurement. The average power measurements were acquired in a slow time
scale, since the power meters have a bandwidth of approximately 50 Hz. For
each $I_{\text{pump}}$, intracavity and extracavity time series, RF spectra
and average powers were acquired. The optical spectrum was measured after the
laser output. All time series were saved directly from the oscilloscope in 200
ns windows and, for each $I_{\text{pump}}$, 20 intracavity and extracavity
windows were acquired and processed. It is important to note that the
photodetectors are much slower than the femtosecond pulses, but the peak
values of the time series are proportional to the pulse energy. It still
faithfully represents the periodic nature of the dynamics as shown in our
previous work [1].

As stated previously, our main interest is in characterizing the nonlinear
transmission through the polarizing beamsplitter. In order to do so, we
measure average powers after its two exit ports, as shown in Fig. 1(b). The
beam reflected by the PBS is the output of the laser and detector D1 measures
its power, as shown in Fig. 1(a). On the other hand, it is not possible
measure the power transmitted through the PBS directly, as this would prevent
the laser from operating. To circumvent this problem we use the zero order
reflection from one of the gratings, collected by detector D4. Careful
measurements of the conversion factors for the detectors were performed and we
determine that the actual powers immediately after the PBS must be multiplied
by 1.3 and 17.9 for D1 and D4, respectively. The internal losses in the PBS
are negligible and the total average power that leaves the gain medium is
given by the from the sum of reflected and transmitted average powers.

##  3 Results and discussion

Using $I_{\text{pump}}$ as a control parameter, different dynamical behaviors
are observed by looking at the variation of the amplitude of the pulses in a
time series. To investigate distinct operating regimes, the evolution of the
RF spectrum was systematically registered for intracavity and extracavity
signals. As both presented similar dynamics, Fig. 2 presents only the
extracavity one. For 525 mA $<I_{\text{pump}}<$ 545 mA, there is only one
single peak at 125 MHz and it corresponds to a regular mode-locked regime,
where pulses have the same amplitude. For 545 mA $<I_{\text{pump}}<$ 829 mA, a
new peak arises at approximately 62.5 MHz, half the fundamental frequency,
corresponding to period doubling. For 829 mA $<I_{\text{pump}}<$ 969 mA,
another peak becomes present at approximately 31.25 MHz, one quarter the
fundamental frequency, corresponding to period quadrupling. All the others
visible frequencies in the RF spectrum are higher order beat notes of these
ones.

Figure 2: Evolution of the extracavity RF spectrum showing period-one, period-
two and period-four. Figure 3: Intracavity (upper row) and extracavity (lower
row) time series (fast time scale measurements): (a) and (e) period-one (529
mA), (b) and (f) period-two (721 mA) and (c) and (g) period-four (957 mA).
Bifurcation diagrams (d) intracavity and (h) extracavity.

The intracavity and extracavity time series also confirm the nonlinear
dynamics observed in the evolution of the RF spectrum. This is seen in the
representative time series shown in Fig. 3, where the upper row refers to the
intracavity signal and the lower row refers to the extracavity one. The
regular mode-locked regime appears in Fig. 3(a) and 3(e), each of them
presenting pulses with constant peak values. However, this uniformity does not
exist either in Fig. 3(b) or 3(f), where the amplitude of the pulses
alternates between two values presenting a period-two regime. Similarly, it
happens to the series presented in Fig. 3(c) and 3(g), where the amplitude is
modulated among four values, showing a period-four regime. Fig. 3(g) draws
attention to the existence of pulses much higher than the others, a very
different pattern from the other series. For every four pulses, the largest of
them has the amplitude almost ten times bigger than the smallest one, making
it more susceptible to the nonlinearity of the gain medium. This pattern
occurs throughout the extracavity period-four regime, which exhibits a much
higher modulation than the intracavity one shown in Fig. 3(c), being relevant
to the analysis because amplitude modulation has a direct relation to the
nonlinear transmission curve. All these regimes are better visualized in the
bifurcation diagrams presented in Fig. 3(d) and 3(h), where the peak values
from the time series are plotted versus $I_{\text{pump}}$ as control
parameter. The branches of the intracavity diagram in Fig. 3(d), remain closer
than the branches of the extracavity in Fig. 3(h). Notice that in Fig. 3(d)
the pulse amplitude seems to saturate, which implies in higher extracavity
amplitudes, as shown in Fig. 3(h). The lack of resolution of the intracavity
bifurcation diagram in Fig. 3(d) is due to the smaller difference between the
intensities of the time series peaks.

Figure 4: (a) Intra, extra and total average power (slow time scale
measurements). (b) Evolution of extracavity optical spectrum. Figure 5:
Nonlinear transmission curve: (a) PBS transmission $T$ versus
$P_{\text{Total}}$, (b) First derivative of $T$ versus $P_{\text{Total}}$.

In order to clarify the meaning of the next set of data, we describe the
relation between the total nonlinear cavity loss, $L(P)$ defined in Eq. (1),
and the nonlinear transmission, $T(P)$ at the PBS. Figure 1(b) shows the
relevant input and output ports of the PBS. In order to achieve mode-locked
operation, the optical feedback must favor higher powers and the portion of
light reflected from the PBS must decrease with power. This means the
nonlinear reflectivity should be of the form $R(P)=R_{\text{PBS}}-\gamma\,P$.
As the internal losses in the PBS are negligible, its nonlinear transmission
is $T(P)=1-R(P)$, which may be rewritten as,

$T(P)=1-R_{\text{PBS}}+\gamma\,P\equiv T_{\text{PBS}}+\gamma\,P\;,$ (2)

which increases with power, as desired. Importantly, we emphasize that both
$R_{\text{PBS}}$ and $\gamma$ may be controlled externally by proper setting
of the waveplates. From Eq. (1), the total intracavity loss in the laser is
given by,

$L(P)=L_{\text{PBS}}+L_{\text{CAV}}-\gamma\;P,$ (3)

where $L_{\text{PBS}}=R_{\text{PBS}}$ and the overall cavity linear loss is
$L_{0}=L_{\text{PBS}}+L_{\text{CAV}}$, where $L_{\text{CAV}}$ represents all
the other linear cavity losses, such as coupling losses back into the fiber,
grating losses, etc.

From Fig. 1(b) the nonlinear reflectivity is obtained from measurements of the
intracavity and extracavity powers, $P_{\text{Intra}}$ and $P_{\text{Extra}}$,
respectively. The total incident power is
$P_{\text{Total}}=P_{\text{Intra}}+P_{\text{Extra}}$ and the nonlinear
reflectivity and transmission are

$\displaystyle R(P)=\frac{P_{\text{Extra}}}{P_{\text{Total}}},$ (4a)
$\displaystyle\quad T(P)=\frac{P_{\text{Intra}}}{P_{\text{Total}}}.$ (4b)

At this point it is important to remind that we use slow detectors for the
power measurements. This means that the power is time-averaged over the pulse
train. In particular, if the laser is in the period-two regime, average values
of the low and high pulse energies are measured and one may ask if the details
of the dynamic regime may be lost. As we will explain, the answer to this
question is no.

We observe a curious phenomenon when the average powers are measured as
function of the pump current: neither $P_{\text{Intra}}$ nor
$P_{\text{Extra}}$ grows linearly with $I_{\text{pump}}$. This is in
contradiction to what would be expected for a usual laser below saturation,
for which the output power $P_{\text{Extra}}$ varies linearly with pump power.
Instead, each curve presents a distinct behavior, especially in the period-
four regime. In Fig. 4(a), for $I_{\text{pump}}<829$ mA the extracavity and
intracavity powers present opposite concavities. When $I_{\text{pump}}$ is
increased above 829 mA, $P_{\text{Intra}}$ remains approximately constant
while $P_{\text{Extra}}$ increases at a higher rate than previously.

A possible explanation is that there is an abrupt variation in the
distribution of instantaneous powers observed in the extracavity time series
where the amplitude of some pulses increases faster with $I_{\text{pump}}$.
The pulses with higher powers access regions of the PBS transmission curve
where the nonlinear effects are more prominent and a description taking into
account contributions beyond the first order power dependence, predicted in
Eq. (2), is necessary [15]. This alters the shape of the effective
transmission, the magnitude of the NPR and the measured average powers. Hence,
the existence of pulses with higher powers in the extracavity signal are
associated with a more pronounced NPR.

By adding $P_{\text{Intra}}$ and $P_{\text{Extra}}$, another interesting
effect becomes evident: the sum shows an approximately linear growth with
$I_{\text{pump}}$, as seen by the curve $P_{\text{Total}}$ in Fig. 4(a). To
appreciate the meaning of this observation we remind that $P_{\text{Total}}$
is the average power immediately after the gain section of the laser, just
before the PBS. If one of the pulses in the pulse train has a higher peak
power, its transmission through the PBS is higher but then the gain it
experiences is smaller, due to gain saturation. The nonlinear contributions to
gain saturation and nonlinear losses seem to compensate each other. The small
signal (unsaturated) gain increases linearly with $I_{\text{pump}}$, and the
net result is that the power after the gain section follows the small signal
gain. In fact, the usual mode-locked operation requires this sort of balance
to establish a stable pulse train. The difference in the present situation is
that it does not occur on each individual pulse, but on the average pulse
power, even during the pulse period doubling sequence.

The transition from period-two to period-four is also followed by changes in
the extracavity optical spectrum, with the emergence of sidebands at 1025 nm,
as shown in Fig. 4(b). These sidebands persist up to $I_{\text{pump}}$ = 969
mA and disappear when $I_{\text{pump}}<$ 829 mA, suggesting that it is
associated with period quadrupling, when pulses with greater intensities
arise. However, for the bifurcation around $I_{\text{pump}}$ = 545 mA, no
transition was observed either in the average power or in the optical
spectrum. This happens because, in the transition from regular regime to
period-two, a small modulation gradually appears in the pulse intensities,
whereas, in the bifurcation to period-four, occurs an abrupt change with
immediate appearance of large pulses, making it possible to access
nonlinearities previously unreachable. Larger intensities are associated with
larger self-phase modulation, which broadens the optical spectrum.

A clear experimental signature of the transitions between different dynamical
regimes is shown in Fig. 5(a), where the nonlinear transmission $T$ as a
function of $P_{\text{Total}}$ is presented. Note that the inclination of the
curve changes around $P_{\text{Total}}$ = 165 mW (or $I_{\text{pump}}$ = 545
mA), which coincides with the period doubling. As pointed out, this
bifurcation was not manifested in individual power measurements or in the
optical spectrum, so that the nonlinear transmission curve brings more
information and appears as a more complete characterization for the study of
nonlinear dynamics of the laser. The transition to period-four occurs at
$P_{\text{Total}}\approx$ 247 mW, which corresponds to $I_{\text{pump}}$ = 829
mA, in agreement with the previous results in Fig. 4. These transitions are
better visualized in the derivative of $T$ with respect to $P_{Total}$ shown
in Fig. 5(b), where abrupt changes occur for each bifurcation, making clear
the existence of three different windows that relate to (I) regular mode-
locked, (II) period-two and (III) period-four regimes. These discontinuities
are a consequence of the changes in the NPR which occur for the different
dynamical regimes.

Many numerical models have been proposed to describe the dynamics of MLFL.
However, accurate quantitative comparison between experiments and theory is
difficult to achieve, due to the complexity of the physical system and the
large number of parameters affecting the resulting dynamics [19, 20]. Any of
those models could be used to show that dynamical systems undergoing
bifurcation display sudden variations in the average value of the dynamical
variables [21]. The average power in lasers is just another example, as seen
for instance in the average power of far infrared gas lasers with saturable
absorber [22]. In this last reference has been shown that the average power
emitted by the laser has a sudden change in its slope as the pump power is
scanned through bifurcations leading to different dynamical regimes. We show a
similar phenomenon in our system in Fig. 5.

##  4 Conclusions

In conclusion, we have measured the internal nonlinear transmission curve of
an Yb MLFL. We have performed measurements while the laser performs
transitions between the standard mode-locking, to period-two and period-four
regime in the pulse energy. It is shown that this curve presents a clear
experimental signature of the dynamical regime transitions where the slope
suddenly varies as the regime changes. We have shown that the evolution of the
intracavity and extracavity average powers are different, due to the variation
of the NPR as the pump power increases. However, $P_{\text{Total}}$ presents a
linear behavior, demonstrating complementary character of $P_{\text{Intra}}$
and $P_{\text{Extra}}$, which is a direct consequence of the energy
conservation. As a next step, we are interested in extending the use of the
nonlinear transmission curve to study harmonic mode-locking and multi-pulsing
regimes and preliminary measurements in this direction have already been
performed.

This work was supported by the Brazilian research agencies CNPq
(441668/2014-3, INCT-IQ 465469/2014-0), CAPES (PRONEX 534/2018,
23038.003382/2018-39), and FACEPE (APQ - 1178 - 1.05/14).

## References

  * [1] Lucas B. A. Mélo and Guillermo F. R. Palacios and Pedro V. Carelli and Lúcio H. Acioli and José R. Rios Leite and Marcio H. G. de Miranda, “Deterministic chaos in an ytterbium-doped mode-locked fiber laser,” Opt. Express 26, 13686-13692 (2018).
  * [2] Antoine F. J. Runge and Claude Aguergaray and Neil G. R. Broderick and Miro Erkintalo, “Raman rogue waves in a partially mode-locked fiber laser,” Opt. Lett. 39, 319-322 (2014).
  * [3] Zhanwei Liu and Shumin Zhang and F. W. Wise, “Rogue waves in a normal-dispersion fiber laser,” Opt. Lett. 40, 1366-1369 (2015).
  * [4] Antoine F. J. Runge and Neil G. R. Broderick and Miro Erkintalo, “Observation of soliton explosions in a passively mode-locked fiber laser,” Optica 2, 36-39 (2015).
  * [5] Grelu, Philippe and Akhmediev, Nail “Dissipative solitons for mode-locked lasers,” Nature Photonics 6 84 (2012).
  * [6] Shian Zhou and Dimitre G. Ouzounov and Frank W. Wise, “Passive harmonic mode-locking of a soliton Yb fiber laser at repetition rates to 1.5 GHz” Opt. Lett. 31, 1041-1043 (2006).
  * [7] Feng Li and P. K. A. Wai and J. Nathan Kutz, “Geometrical description of the onset of multi-pulsing in mode-locked laser cavities,” J. Opt. Soc. Am. B 27, 2068-2077 (2010).
  * [8] Komarov, Andrey and Leblond, Hervé and Sanchez, François, “Multistability and hysteresis phenomena in passively mode-locked fiber lasers,” Phys. Rev. A 71, 053809 (2005).
  * [9] Haboucha, Adil and Komarov, Andrey and Leblond, Hervé and Sanchez, François and Martel, Gilles, “Mechanism of multiple pulse formation in the normal dispersion regime of passively mode-locked fiber ring lasers,” Optical Fiber Technology 14, 262-267 (2008).
  * [10] Andy Chong and William H. Renninger and Frank W. Wise, “All-normal-dispersion femtosecond fiber laser with pulse energy above 20nJ,” Opt. Lett. 32, 2408-2410 (2007).
  * [11] , Bülend Ortaç and Martin Baumgartl and Jens Limpert and Andreas Tünnermann, “Approaching microjoule-level pulse energy with mode-locked femtosecond fiber lasers,” Opt. Lett. 34, 1585-1587 (2009).
  * [12] Nelson, L.E. and Jones, D.J. and Tamura, K. and Haus, H.A. and Ippen, E.P., “Ultrashort-pulse fiber ring lasers,” Applied Physics B 65, 277-294 (1997).
  * [13] H. A. Haus and J. G. Fujimoto and E. P. Ippen, “Structures for additive pulse mode locking,” J. Opt. Soc. Am. B 8, 2068-2076 (1991).
  * [14] M. Hofer and M. E. Fermann and F. Haberl and M. H. Ober and A. J. Schmidt, “Mode locking with cross-phase and self-phase modulation,” Opt. Lett. 16, 502-504 (1991).
  * [15] N. R. Newbury and B. R. Washburn, “Theory of the frequency comb output from a femtosecond fiber laser,” IEEE Journal of Quantum Electronics 41, 1388-1402 (2005).
  * [16] Christopher R. Jones and J. Nathan Kutz, “Stability of mode-locked pulse solutions subject to saturable gain: computing linear stability with the Floquet–Fourier–Hill method,” J. Opt. Soc. Am. B 27, 1184-1194 (2010).
  * [17] Akhmediev, N. and Soto-Crespo, J. M. and Town, G., “Pulsating solitons, chaotic solitons, period doubling, and pulse coexistence in mode-locked lasers: Complex Ginzburg-Landau equation approach,” Phys. Rev. E 63, 056602 (2001).
  * [18] Wei, H. and Li, B. and Shi, W. and Zhu, X. and Norwood, R. A. and Peyghambarian, N. and Jian, S. “General description and understanding of the nonlinear dynamics of mode-locked fiber lasers,” Scientific Reports 7, 1292 (2017).
  * [19] Shaokang Wang and Andrew Docherty and Brian S. Marks and Curtis R. Menyuk, “Boundary tracking algorithms for determining the stability of mode-locked pulses,” J. Opt. Soc. Am. B 31, 2914-2930 (2014).
  * [20] Shaokang Wang and Brian S. Marks and Curtis R. Menyuk, “Comparison of models of fast saturable absorption in passively modelocked lasers,” Opt. Express 24, 20228-20244 (2016).
  * [21] Cavalcante, H. L. D. S. and Rios Leite, J. R., “Bifurcations and averages in the logistic map,” Dynamics and Stability of Systems 15, 35-41 (2000)
  * [22] Oliveira-Neto, L. B. and da Silva, Glauber J. F. T. and Khoury, A. Z. and Rios Leite, J. R., “Average intensity and bifurcations in a pulsed laser,” Phys. Rev. A 54, 3405-3408 (1996).