CoOMBE: A suite of open-source programs for the integration of the optical Bloch equations and Maxwell-Bloch equations

The programs described in this article and distributed with it aim (1) at integrating the optical Bloch equations governing the time evolution of the density matrix representing the quantum state of an atomic system driven by laser or microwave fields, and (2) at integrating the 1D Maxwell-Bloch equations for one or two laser fields co-propagating in an atomic vapour. The rotating wave approximation is assumed. These programs can also be used for more general quantum dynamical systems governed by the Lindblad master equation. They are written in Fortran 90; however, their use does not require any knowledge of Fortran programming. Methods for solving the optical Bloch equations in the rate equations limit, for calculating the steady-state density matrix and for formulating the optical Bloch equations in the weak probe approximation are also described.


PROGRAM SUMMARY
Program Title: CoOMBE CPC Library link to program files: (to be added by Technical Editor) Developers' repository link: https://github.com/durhamqlm/CoOMBECode Ocean capsule: (to be added by Technical Editor) Licensing provisions: GPLv3 Programming language: Fortran 90 Nature of problem: The present programs can be used for the following operations: (1) Integrating the optical-Bloch equations within the rotating wave approximation for a multi-state atomic system.At the choice of the user, the calculation will return either the time-dependent density matrix at given times or the density matrix in the long time limit if the system evolves into a steady state in that limit.The calculation can be done with or without averaging over the thermal velocity distribution of the atoms.The number of atomic states which can be included in the calculation is limited only by the CPU time available and possibly by memory requirements.An arbitrarily large number of laser or microwave fields can be included in the calculation if these fields are all CW.This number is currently limited to one or two for fields that are not all CW.The calculation can be done in the weak probe approximation, or in the rate equations approximation, or without assuming either of these two approximations.Calculating refractive indexes, absorption coefficients and complex susceptibilities is also possible.(2) Integrating the 1D Maxwell-Bloch equations in the slowly varying envelope approximation for one or two fields co-propagating in a single-species atomic vapour.Although geared towards the case of atoms interacting with laser fields, this code can also be used for more general quantum systems with similar equations of motion (e.g., molecular systems, spin systems, etc.).Solution method: The Lindblad master equation is expressed as a system of homogeneous first order linear differential equations,

Introduction
The programs described in this article have been developed for modelling general atomic or molecular systems interacting with one or several laser or microwave fields resonant or nearly resonant with atomic transitions, the interaction being treated within the rotating wave approximation.They can also be used for more general quantum systems, e.g., spin systems, governed by similar equations of motion.Their main focus is on the calculation of the density matrix representing the state of the system as obtained by integrating the optical Bloch equations (i.e., the Lindblad master equation for such systems).The populations and coherences can be calculated as time-dependent functions.Alternatively, for systems driven by CW fields and evolving to a steady state, they can also be obtained in the long time limit.Calculations of complex susceptibilities, refractive indexes and absorption coefficients are also possible.The present programs were originally developed for studying the formation of optical solitons and two-colour quasisimultons in an optical vapour, which required the integration of the 1D Maxwell-Bloch equations in the slowly envelope approximation [1].The Maxwell-Bloch solvers developed at that occasion are included in this library as they are of general interest and share the same user interface. 1 A number of open source programs are already available for tackling similar or related calculations, namely general purpose programs for the modelling of open quantum systems, programs more specific to Atomic Physics calculations, and programs solving the Maxwell-Bloch equations in various approximations.The general purpose programs include, in particular, Qutip [2,3], written in Python, its predecessor, Quantum Optics Toolbox [4,5], written in MATLAB, and a more recent alternative, QuantumOptics, written in Julia [6].They also include several MATLAB programs primarily intended for educational purposes [7,8,9], two Quantum Monte Carlo programs written in C++ [10,11] and a MATLAB program focusing on the optimal control of the dynamics of quantum systems interacting with external electromagnetic fields [12].General programs solving the optical Bloch equations for atomic systems have also been published, including the Atomic Density Matrix package [13], which is written in Mathematica and also supports more general quantum optical calculations, a collection of Python tools for modelling few-level atom-light interactions [14], and PyLCP [15], a Python program oriented towards the modelling of laser cooling but also allowing for general solutions of the optical Bloch equations.The Elecsus program, also written in Python, is specialised to the case of an atomic vapour addressed by a single probe field and offers powerful facilities for the analysis of experimental absorption spectra for that particular case [16,17].The necessary atomic data are provided by Atomic Density Matrix, PyLCP and Elecsus for species of current interest.The Maxwell-Bloch solvers include mbsolve [18], a C++ program for the 1D propagation of a field in the plane wave approximation or a field confined to a wave guide, and QuEST [19,20,21,22], also written in C++, which was developed for modelling the interaction of an electromagnetic field with multiple 2-level quantum dots in 3D.Neither mbsolve nor QuEST assume the slowly varying envelope approximation; mbsolve does not assume the rotating wave approximation either and is not restricted to a 2-state medium (the Maxwell-Bloch solver included in the present package is simpler: it takes the medium to be homogeneous and assumes both the slowly varying envelope approximation and the rotating wave approximation, which is appropriate, e.g., in calculations of self-induced transparency for many-cycle pulses).
The present programs may nonetheless be of interest in view of their generality, their scalability to large systems, and the wide choice of integration methods they offer.While written in Fortran 90 for speed and convenience, no knowledge of this language is necessary for using them.It is expected that they will be further extended in the future, in particular by coupling them with programs providing the atomic data required for calculations on atomic systems of current experimental interest.Future versions will be published at the URL https://github.com/durham-qlm/CoOMBE,a GitHub repository 1 Python codes developed in the course of the work reported in [1] are published at the URL https://github.com/tpogden/maxwellbloch. of the Quantum, Light and Matter research group of Durham University.
General information about the distribution is given in Section 2. The computational methods implemented in these program and the theoretical framework are outlined in Section 3. Information about installing and using this software and information can be found in Section 4 and also, at much greater lengths, in the User Manual.The reuse of codes written by other authors in the present programs is acknowledged in Section 5.The main text is accompanied by a number of more technical appendices: the reduction of Lindblad master equation to rate equations is explained Appendix A, the calculation of a steady-state density matrix in the long time limit in Appendix B, and the implementation of the weak probe approximation in the present computational framework in Appendix C. Useful results concerning weak probe calculations for a single field are gathered in Appendix D. Examples of the use of these programs are given in Appendix E and Appendix F, and further examples can be found in the examples directory included in the distribution.Advice about how to run this software without the installation of a Fortran compiler and supporting libraries, through a Podman container [23], are provided in Appendix G and in the GitHub repository.
We are not aware that the methods described in Appendix A, Appendix B and Appendix C are widely known or previously published.

Organisation into program units
This library contains several modules and external subroutines, as follows: • The general settings modules.This module sets several key parameters, in particular the variable nst which defines the number of states to be considered in the calculation.More information about this module and these parameters can be found in Section 4.2.
• The obe constants module, which defines fundamental physical constants used elsewhere in the code [24].
• The obe module, which forms the main part of the library.It contains a number of subprograms, many of which are private to this module (i.e., cannot be called from outside the module).These subprograms are concerned with solving the optical Bloch equations and/or forming a userfriendly interface with the ldbl module.
• The mbe module, grouping program units concerned with solving the Maxwell-Bloch equations and with solving the optical Bloch equations for time-dependent fields.
• The ldbl module, which contains a number of subprograms concerned with setting up and solving the Lindblad master equation.This module is the core of the library.However, the subprograms it contains can be accessed more conveniently through subprograms forming part of the obe or mbe modules, and using those does not require any knowledge of the inner working of the ldbl module.For this reason, the latter is not addressed in the present article.The reader is referred to the detailed documentation for general information about its content.
• The external subroutine ext setsys, which is used only for communicating information between the obe and mbe modules.
• The external subroutines fcn dummy and solout dummy, which are provided for compatibility with the original code of the DOP853 ODE solver mentioned in Section 3.1.3.
• The ldblstore module, which is used to store certain intermediate results produced by programs contained in the ldbl module.
• The driveall program, described in Section 4.4, which offers a simple interface with obe and mbe and makes it possible to use these codes without any Fortran programming.
These various components are grouped into files as outlined in Table 1.
Besides a number of program units intended for internal use only, the obe and mbe modules currently contain a total of 35 user-facing subprograms (Table 2), i.e., subprograms providing an interface between the internal program units and a userwritten driving program.

Documentation and examples
The distribution includes a User Manual complementing the present article.This document contains further information about all the programs units forming this library, a detailed description of the user-facing subprograms, detailed information about the use of the driveall program, and a short tutorial explaining how the Hamiltonian of an atomic system interacting with an electromagnetic field treated in the rotating wave approximation can be cast into the form of Eq. (11).
Further examples illustrating the use of various features of this library are also provided.The corresponding files and documentation can be found in the examples directory.

General formulation
The obe and mbe codes have been developed for modelling atomic systems driven by laser fields or other coherent electromagnetic fields and composed of two or more atomic states, the fields being resonant or close to resonance with transitions between these states.The codes can also be used to calculate the density matrix for more general N-state quantum systems interacting with a superposition of M electromagnetic fields, as long as the rotating wave approximation can be assumed.The general settings module.

ldbl.f90
The ldbl and ldblstore modules and the fcn dummy and solout dummy subroutines.

obe.f90
The obe and obe constants modules and the ext setsys subroutine.

driveall.f90
The driveall program.examples Examples of the use of this software, including the files and the program listed in Appendix E and Appendix F.
Each field is described by a real electric field vector, E α (r, t), the total electric field of the applied light at position r and time t being E(r, t), with The calculation assumes that each of the E α (r, t)'s can be written as the product of a slowly-varying envelope and a planewave carrier.Specifically, where k α is the wave vector of field α and ω α is its angular frequency (ω α > 0).The field amplitudes E α may be complex and may vary in time.The polarisation vectors are assumed to be constant and of unit norm: Given Eqs. ( 2) and (3), the intensity of a continuous wave (CW) field is related to its complex amplitude by the following equation: This relation generalizes to the case of a pulsed field of envelope E(t), provided the pulse encompasses more than a few optical cycles: ϵ 0 c |E(t)| 2 /2 can be taken to be the instantaneous intensity at time t.Intensities are easily converted into electric field amplitudes and conversely by making use of the fact that an intensity of exactly 1 mW cm −2 corresponds to an electric field amplitude of 86.8021 V m −1 .
The states coupled to each other by the field(s) are assumed to be orthonormal eigenstates of the field-free Hamiltonian, Ĥ0 .We denote these states by the ket vectors |i⟩, i = 1, . . ., N, and the corresponding field-free eigenenergies by ℏω (i) : Typically, these N states form two or more groups differing considerably in energy and each of the fields is resonant or close to Specifies convergence criteria for the DOP853 ODE solver.

Auxiliary routines obe coher index
Given the indexes of two different atomic states, i and j, returns the indexes of the components corresponding to the real and imaginary parts of the coherence ρ i j in the 1D representation of the density matrix.obe fieldtocfield Given a variable of type obefield containing the details of a field, returns a variable of type obecfield containing the same details.obe find campl Given a complex dipole moment and a complex Rabi frequency, returns the corresponding complex electric field amplitude.obe find rabif Given a complex dipole moment and a complex electric field amplitude, returns the corresponding complex Rabi frequency.

obe init rho
Returns the density matrix of a mixed state with given populations and zero coherences.

obe pop index
Returns the index of the component corresponding to a specified population in the 1D representation of the density matrix.

obe susceptibility
Given the relevant coherences, calculates the complex susceptibility, refractive index and absorption coefficient.
resonance with transitions between states of one of these groups and states of one of the other groups.E.g., in rubidium, these groups could be the 5S 1/2 (F, m F ) states, the 5P 1/2 (F, m F ) states, the 5P 3/2 (F, m F ) states, etc., and the calculation could involve a field resonant or nearly resonant on a transition between one of the 5S 1/2 (F, m F ) states and one of the 5P 1/2 (F, m F ) states.These groups of energetically close states are denoted by G 1 , G 2 , . . ., G K in the following.The states belonging to a same group may or may not differ in energy, depending on the system.Either way, the energies ℏω (i) of the states belonging to a same group can be referred to a reference energy, ℏω ref , from which each one differs by an energy offset ℏδω (i) : for group k, A reference energy ℏω ref (k) could be, for example, the energy of one of the basis states, or the centroid of a group of hyperfine levels.More generally, ℏω ref (k) can be any energy reference appropriate for the problem at hand.The calculation also assumes that the interaction with the fields is taken into account within the electric dipole approximation, which amounts to neglecting the spatial variation of E(r, t).For simplicity, the vector r is taken to be zero in the following.The exp(±i k α • r) phase factors can be subsumed into the complex amplitudes E α should they be relevant.The Hamiltonian of the system thus takes on the following form: where D is the atom's dipole operator.In terms of the relevant position operator, X, where e is the absolute charge of the electron (e > 0).The matrix elements of the operator εα • X would typically be obtained as the product of a reduced matrix element and an angular factor.The corresponding complex Rabi frequencies Ω α;i j are defined as follows throughout the code: with the convention that Ω α;i j = 0 if states i and j are deemed not to be coupled by field α, e.g., because this transition would be excessively far from resonance.It should be noted that this definition of the complex Rabi frequency includes the (negative) −e factor multiplying the position operator.It may differ, in sign and otherwise, from the definition of the Rabi frequency used by other authors.This Hamiltonian cannot be treated in its full complexity by the present software.Rather, the obe and mbe routines are based on a simplified Hamiltonian, Ĥ′ , derived from Ĥ(t) by neglecting any excessively far detuned transition, making the rotating wave approximation and passing to slowly varying variables by a unitary transformation.This transformed Hamiltonian is assumed to have the following general form: 11) where ∆ α is the frequency detuning of field α, the a iα 's are numerical factors, and What the factors a iα are and how the frequency detunings are defined in terms of the energies of the relevant states and the angular frequencies ω α varies from system to system, as is explained in Appendix A of the User Manual.For instance, for the 3-state system considered in Appendix E, and, as can be seen from Eq. (E.1), a 11 = a 12 = a 22 = 0 and a 21 = a 31 = a 32 = −1.We stress that these frequency detunings, as defined, are angular frequencies, like the Rabi frequencies Ω i j .It can be noted that Ĥ′ is a self-adjoint operator, as expected, since Ω ji = Ω * i j within the above definition of the Rabi frequencies.For most systems, Ĥ′ is constant in time if all the fields considered are CW fields.
The optical Bloch equations are the equations of motion for the individual components of the density matrix for an open quantum system interacting with classical electromagnetic fields.They are obtained from the Lindblad master equation, where ρ is the density operator describing the state of the system and the Ĉn 's are certain operators called jump (or collapse) operators.The latter include the operator Γ i j | i ⟩⟨ j | if state j relaxes to state i at a rate Γ i j by spontaneous decay or some other mechanism.It is customary to add phenomenological terms in decay at an additional rate γ i j due to pure dephasing effects such as collisional broadening.
The obe and ldbl modules calculate the density matrix, ρ, representing the density operator ρ in the {|i⟩} basis -i.e., the elements of ρ are the matrix elements of ρ: Re ρ i j = Re ρ ji and Im ρ i j = −Im ρ ji since the density matrix is Hermitian.These relations are used within the obe, mbe and ldbl modules to store and calculate this matrix as a column vector of N 2 real numbers, r, rather than as a 2D array of N 2 complex numbers.Specifically, if the states are labelled 1, 2, 3,. . .as done throughout this section, Accordingly, the Lindblad equation is recast as a set of homogeneous linear relations between the elements of r and the elements of ṙ, the time derivative of r: where L is a N 2 × N 2 real matrix.Much of the obe and ldbl code aim at constructing this matrix given the parameters of the system and at integrating Eq. ( 18), either as written or after further transformation.(Readers interested in knowing the details of how the matrix L is constructed are referred to the information given in the code of the subroutine ldbl reformat rhs cmat contained in the subroutine ldbl set rhsmat of the ldbl module.)

Inhomogeneous broadening
The obe and mbe modules make it possible to take inhomogeneous broadening into account in the calculation of the density matrix.The codes are specifically geared towards the case of Doppler broadening arising from the free thermal motion of atoms in an atomic vapour.However, they can be easily generalised to other cases of Gaussian broadening if required.Extending them to Doppler broadening for a non-Maxwellian distribution of atomic velocities is also possible.
In the current state of development of the obe and mbe modules, Doppler averaging is possible only for co-propagating or counter-propagating co-linear fields.The internal state of an atom depends on the component of its velocity vector in the direction of propagation of the field, v, owing to the Doppler shift of the detunings ∆ α .To first order in 1/c, if the wave vector k α is oriented in the positive z-direction or if it is oriented in the negative z-direction.Correspondingly, the matrix L appearing in Eq. ( 18) depends on v, and so does the solution vector r.Averaging the latter over the Maxwellian distribution of atomic velocities gives the Doppler-averaged density matrix, ρ av , here represented by the column vector r av : with In this last equation, u is the rms velocity of the atoms in the z-direction: u = √ 2k B T/M, where k B is Boltzmann constant, T is the temperature of the vapour and M is the mass of the atom.
The obe and mbe modules include code calculating the integral over v either by numerical quadrature or by expressing the integral in terms of the Faddeeva (or Faddeyeva) function, w(z) [25]: In terms of the complementary error function [25], The approach based on the Faddeeva function applies only to Doppler averaging of the steady state density matrix.It is outlined in Appendix B and Appendix D.
The numerical quadrature method is more general.The quadrature abscissas {v k } and quadrature weights {w k } used by the obe and mbe modules can either be provided by the user or calculated internally.As the program sets the quadrature weights should not include the velocity distribution f M (v).Since the Doppler effect is taken into account to first order in 1/c only, as per Eqs.( 19) and ( 20), the matrix L varies linearly with v: where L 0 and L 1 do not depend on v.These two matrices are easily constructed, which makes Eq. ( 26) an efficient way of recalculating L for each value of v. Replacing f M (v) by another velocity distribution, should this be necessary, would only require minor changes to the codes.

Integrating the optical Bloch equations
Integrating Eq. ( 18) subject to specified initial conditions gives the density matrix as a function of time.Unless the size of the system is excessively large, this operation is amenable to standard numerical methods.This library provides five subroutines to this effect, namely obe Doppler av td A, obe Doppler av td B and obe tdint, for CW fields, and mbe tdint 1 and mbe tdint 2, for fields with a timedependent complex amplitude E α (t).The obe routines can handle an arbitrary number of applied fields, whereas mbe tdint 1 and mbe tdint 2 are respectively limited to one and two fields.Both obe Doppler av td A and obe Doppler av td B calculate the Doppler-averaged time-dependent density matrix.These two routines differ by their memory and CPU times requirements.The obe tdint routine calculates the timedependent density matrix without Doppler averaging.Doppler averaging is optional for the two mbe routines.
Each of these five routines offers a choice of integrator between the classic fourth-order Runge-Kutta method, Butcher's fifth-order Runge-Kutta method [26] and an adaptive ODE integrator (the DOP853 routine of Hairer et al, which is a Dormand-Prince implementation of an explicit eighth-order Runge-Kutta method [27,28]).A solution based on the right and left eigenvectors of the matrix L is also implemented, and can be contemplated if this matrix is time-independent (which is normally the case if the applied fields are CW).These eigenvectors fulfill the equations and with N = N 2 .In many cases of interest, the initial density matrix vector can be written as a linear combination of the v j 's.I.e., there exist complex coefficients c 1 , c 2 , . . ., c N such that In this case, the density matrix can be obtained for all times as However, the existence of such a set of coefficients is not guaranteed since the matrix L is not symmetric and may be defective.The subprogram obe tdint offers the option to attempt to expand r(t = t 0 ) as per Eq. ( 29) with and if this attempt is successful (it normally is), use Eq. ( 30) to propagate the density matrix in time.

Rate equations
The optical Bloch equations can be transformed into a smaller system of rate equations if the elements of the density matrix can be divided into two classes, R and S, depending on whether they converge to steady values much more rapidly (R) or much more slowly (S) than the elements belonging to the other class.Class R typically includes most or all the coherences, class S the populations and, if any, the coherences not included in R.This dichotomy makes it possible to reduce the number of coupled differential equations by adiabatic elimination of the elements belonging to R. The details of this approach can be found in Appendix A.
The routines obe Doppler av td A and obe tdint can solve Eq. ( 18) within this approximation for a superposition of CW fields, with or without Doppler averaging.As the code is written, the set S of the elements of ρ which are actually propagated in time includes all the populations and none of the coherences.The latter are derived from the former through Eq. (A.5) of Appendix A. Time propagation thus involves solving a system of only N coupled differential equations, which is a considerable reduction from the original system of N 2 equations.

Steady state solutions
In many cases, but not all cases, the populations and coherences settle to constant values as time increases if the fields are CW.Then r → r st for t → ∞, where ṙst = L r st = 0. ( The steady-state density matrix represented by the column vector r st is thus an eigenvector of the matrix L corresponding to a zero eigenvalue, and can usually be calculated as such.The calculation follows the same lines as the calculation of r(t) by the eigenvalue method described in Section 3.1.3,except that here only the eigenvectors v j belonging to a zero eigenvalue are included in Eq. ( 30).The optical Bloch equations have no steady state solution if some of the eigenvalues λ j are imaginary.The obe module also supports a different way of obtaining the steady-state density matrix, which is based on transforming the eigenvalue equation L r = 0 into an inhomogeneous system of linear equations, where L ′ is a (N − 1) × (N − 1) square matrix and b is a (N − 1)component column vector.The matrix L ′ and the column vector b are derived from L by a straightforward rearrangement process.The transformation is normally possible due to the unit trace property of the density matrix, which constraints the solutions of this eigenvalue equation.The vector r st representing the steady state density matrix is identical to the solution vector r ′ , apart from one population which can be calculated readily as a linear combination of the other populations.The reader is referred to Appendix B for the details of the method.Calculating r st in this way may be faster than by using the eigenvalue method but will fail if Eq. ( 32) has more than one solutions.It would then be necessary to specify the density matrix that r st develops from in order to obtain a unique solution, which is not overly difficult in the eigenvalue method -and is implemented in the obe module -but would considerably complicate the calculation based on the linear equations method.
Finding the steady state as per Eq. ( 33) also makes it possible to Doppler average the density matrix semi-analytically, as an alternative on the entirely numerical approach mentioned in Section 3.1.2.This semi-analytical route may lead to substantial savings in CPU time as compared to a numerical quadrature.Its principles are outlined in Appendix B.
As pointed out in that appendix, significant savings may also be achieved, along similar lines, in computations involving the calculation of the density matrix for multiple values of a same detuning.
Calculations of the steady-state density matrix are possible only for CW fields.Several routines are provided to this end, namely, for general systems, obe steadystate (for calculations without Doppler averaging), obe Doppler av st (for calculations with Doppler averaging performed semi-analytically as described in Appendix B, and obe Doppler av st numerical (for calculations with Doppler averaging performed by a numerical quadrature).The library also includes a routine specialised to the case of 2-state systems driven by a single field (obe 2state), one specialised to the case of multi-state systems driven by a single field with the calculation organised as explained in Appendix B (obe steadystate onefld), and, as described in next section, several routines specialised to calculations in the weak probe approximation.obe steadystate and obe Doppler av st numerical can handle calculations using the eigenvalue method, which makes it possible to address cases for which the steady state depends on the initial populations.The subroutine obe Doppler av st numerical is normally less efficient than obe Doppler av st.

The weak probe approximation
These programs offer the option of solving the optical Bloch equations within the approximation where one of the fields is considered to be too weak to cause any appreciable optical pumping over the relevant time scales.A calculation within this approximation amounts to calculating the density matrix to first order in the weak field and to all orders in any of the other fields in the problem.The populations are not affected by the former in this case, while the coherences depend linearly on its amplitude, without any power broadening.
The weak field is referred to as the probe field in many applications of these methods, and the weak field approximation as the weak probe approximation.This terminology is also used here.How this approximation is implemented within the obe module is explained in Appendix C.
The calculation of the steady state for a ladder system by the linear equations method may be problematic in the weak probe approximation.Ladder systems here refer to systems in which a set of low energy states, which are the only ones initially populated, are coupled to states of higher energy only by the probe field.The populations of the lower energy states do not vary in time in the weak probe approximation for such systems, and the populations of the higher energy states remain identically zero at all times.The steady state populations are thus the same as the initial ones, which are specified by the user.It is therefore possible to find the steady state coherences by an application of the rate equations method.Referring to Appendix A, the calculation simply amounts to solving Eq. (A.5) for the vector r R , with R including all the coherences and S all the populations.Within this approach, folding the result on a Maxwellian distribution of atomic velocities can also be done in terms to the Faddeeva function, following the same method as outlined in Appendix B but here starting from Eq. (A.5) rather than from Eq. (33).
Calculations within the weak probe approximation can normally be handled by the general computational routines contained in the obe module.However, steady state calculations for ladder system are best done by the subroutine obe steadystate ladder.
The steady state density matrix takes on a particularly simple form in the weak probe approximation if the system comprises only two states or two groups of states coupled by a single field.This case is described in Appendix D. Three specialised routines are provided for tackling such systems, namely obe steadystate onefld weakprb, obe weakfield (a stand-alone subprogram which also calculates the complex susceptibility, the refractive index and the absorption coefficient), and obe steadystate onefld powerbr (for systems in which power broadening is significant but optical pumping is not, as explained in the detailed description of this subprogram given in the User Manual).

The complex susceptibility
Let P(t) be the polarisation generated in the medium by the optical field described by Eq. ( 1).(As mentioned above, we set r = 0 in this equation.The exp(±i k α • r) phase factors are assumed to be taken into account through the complex amplitudes E α if they are relevant.)In terms of a complex susceptibility χ(ω α ), where the . . .stand for contributions oscillating at frequencies other than ω α , if any is present.It is assumed, in the following, that these additional contributions are negligible or absent.
For such systems, where N d is the number density and, as in Section 3.1.1,D is the dipole operator.For each field, the summation runs over all the states |i⟩ and all the states | j⟩ dipole-coupled to each other by this field and such that ℏω (i) > ℏω ( j) with ω (i) − ω ( j) ≈ ω α .This equation can also be written in the following form, which is the one implemented in the programs: The coherences ρ i j 's and therefore the susceptibility χ(ω α ) generally depend on the intensity of all the fields included in the calculation -with the important exception of systems containing only one field and this field is treated within the weak probe approximation (see Appendix D).
Besides the complex susceptibility, the programs can also calculate the corresponding refractive index, n(ω α ), and absorption coefficient, α(ω α ).Here [32], The library contains two routines calculating these quantities, namely obe susceptibility, which uses pre-calculated coherences, and obe weakfield, which is self-contained and computes the necessary coherences within the weak probe approximation for multi-state single-field systems.

The Maxwell-Bloch equations 3.2.1. General formulation
The mbe module addresses the case of a single field or a superposition of two different fields, i.e., a probe field and a coupling field, (co)propagating in the positive z-direction.Solving the Maxwell-Bloch equations for more than two fields or in another geometry is not yet supported.
In general, the spatial and temporal variation of the electric field component of the electromagnetic wave is governed by the equation where P is the medium polarisation and µ 0 is the vacuum permeability.The plane wave approximation is assumed in the calculation performed by the mbe codes.I.e., it is assumed that E and P are constant in any plane perpendicular to the z-axis.
These fields thus depend only on z and t, and the 3D wave equation reduces to the 1D equation This equation can be simplified further, to by making the ansatz and taking into account that the complex amplitudes E α (z, t) and P α (z, t) vary slowly compared to the carriers.As noted above, the library only supports calculations for a single field (M = 1) or a superposition of two fields (M = 2).The field with α = 1 is referred to as the probe field and the field with α = 2 (if present) as the coupling field.
The relationship between the medium polarisation and the state of the atoms is considered in Section 3.1.7,from which it follows that where N d is the medium number density and the summation runs as in Eqs. ( 35) and (36).Changing the time variable t to the shifted time t ′ , with further simplifies Eq. ( 41) to where E α and the coherences ρ i j are now functions of z and t ′ rather than functions of z and t.This last equation governs the propagation of the fields through the medium, as calculated by the present programs.

Implementation
The subroutines mbe propagate 1 and mbe propagate 2 solve Eq. ( 46), with the coherences obtained by solving Eq. ( 15), respectively for the case of a single field (α = 1) or a superposition of two fields (α = 1, 2).The calculation yields the density matrix describing the state of the medium, ρ(z, t ′ ), and the complex amplitude(s) of the propagated field(s), E α (z, t ′ ).These results are calculated on a two-dimensional mesh of values of z and t ′ .The grid points in the z-direction extend from z = z 0 = 0 (the entrance of the medium) to z = z max and are separated by a constant step h: z = z i = z 0 + ih, i = 0, . . ., N z , with z 0 = 0 and h = z max /N z .The distance z max and the number of spatial steps, N z , are set by the user.The complex amplitude of the fields at z 0 must be provided on a mesh of N t + 1 values of t, namely at t = t k with k = 0, 1,. . ., N t .The same mesh is used by mbe propagate 1 and mbe propagate 2 for the shifted time t ′ .Namely, at all z, the grid points in the t ′direction are taken to be at The calculation alternates at each spatial step between obtaining the coherences ρ i j (z, t ′ ) given the field(s) and propagating the field(s) to the next step given these coherences.If Doppler averaging is required, the coherences are obtained for a number of velocity classes and their average, weighted by the Maxwellian velocity distribution, is calculated by numerical quadrature.
The density matrix is calculated at each z i by integrating the optical Bloch equations, starting, at t ′ = t ′ 0 , with initial values determined by the user.A fourth order Runge Kutta rule is used to this end for the integration in time and a predictor-corrector method combining the third order Adams-Bashford rule and the fourth order Adams-Moulton rule for the integration in space.Other choices of methods are also offered.More information about the different possibilities can be found in the User Manual.

Installation
The most recent version of these modules can be found at the URL https://github.com/durham-qlm/CoOMBE.Installing this software only involves downloading the general settings.f90,obe.f90, mbe.f90, ldbl.f90 and driveall.f90files, and editing the general settings.f90file as required (see Section 4.2).The latter is the only program file which may need customisation.
A calculation using these modules requires a driving program, which could be either the driveall program provided in the driveall.f90file or a user-written bespoke Fortran 90 program.The driveall program is described in Section 4.4.Information relevant for the development of a bespoke driver can be found in Section 4.5.The User Manual included in this distribution contains detailed information about the use of driveall and (for Fortran programmers) the use of the various user-facing subroutines contained in these modules.
Running this software requires compiling the programs and linking it to the LAPACK [33] and BLAS [34] libraries.If a Fortran compiler and these two libraries are already installed, compiling these programs could be done, e.g., by the command2 gfortran general_settings.f90ldbl.f90obe.f90 mbe.f90 driveall.f90-llapack -lblas and similarly for a bespoke program.In the latter case, the mbe.f90 file does not need to be compiled if the program does not call any of the mbe subroutine listed in Table 2.
We also provide advice, in Appendix G, for compiling and running these programs through a container, specifically a Podman container [23].This alternative, self-contained method allows the software to be used without installing a Fortran compiler or any supporting libraries directly to the user's machine.This feature is offered for convenience to users not familiar with compiling Fortran codes and does not limit the scope of the program or its output.The distribution includes the files required for running all the examples provided in the examples directory in this way.

Key parameters
The following parameters are defined in the general settings module and must be adapted to the requirements of the intended calculation before compilation: nst: An integer constant which must be given a value equal to the number of states in the model, N. Changing the value of nst is the only editing which may be required across all the modules for adapting the Fortran code to the problem at hands.
kd: An integer defining the kind of many of the variables used in obe and mbe -i.e., defining whether these variables are of real or double precision type (complex or double complex for variables storing complex numbers).Selecting a kind parameter corresponding to real variables rather than to double precision variables will reduce memory requirements and computation time but may also result in larger numerical inaccuracies.
nmn: An integer constant defining how the states are numbered by the user, as explained in Section 4.6.

Required data
The routines provided require various input data, which are problem-dependent and need to be prepared separately.These will typically include: • the energy offset ℏδω (i) defined by Eq. ( 7) for each of the states considered; • the rates of spontaneous decay Γ i j from a state j to a state i, for all the states considered; • any additional dephasing rate γ i j that would need to be included in the calculation to take into account the frequency widths of the fields and/or other pure dephasing effects; • the detuning ∆ α , complex field amplitude E α and wavelength of each of the fields considered, as well as the transition dipole moments for each of the transitions driven by these fields or the corresponding Rabi frequencies; • the initial populations (i.e., the initial values of the diagonal elements of the density matrix); • the temporal profile of the applied field(s), unless these fields are CW or their profile can be calculated internally; • the frequency widths of the fields, if these widths should be taken into account otherwise than through the rates γ i j ; The starting number in the indexing of the states.

filename controlparams
The name of the controlparams file.

Optional parameters filename defaultdata
The name of the defaultdata file.icmplxfld Parameter indicating whether the field amplitudes, dipole moments and Rabi frequencies are specified as real numbers or as complex numbers in the input files.
• the atomic number and the wavelength of each field density in the case of a propagation calculation or a calculation of the complex susceptibility.
All energies and angular frequencies are to be provided as frequencies specified in MHz.E.g., the energy offset ℏδω (i) needs to be provided as the frequency δω (i) /(2π).Wavelengths are to be expressed in nm, densities in number of atoms per m 3 , dipole moments in C m and electric field amplitudes in V m −1 .
Besides the wavelength of each field, calculations involving Doppler averaging will also require the r.m.s.thermal speed of the atoms in the laser propagation direction, u, in m s −1 , and the abscissas and integration weights for the numerical quadrature over the atoms' velocity distribution (unless the calculation uses one of the quadrature rules provided by obe or the integration is done analytically using the Faddeeva function).

Running these codes through the driveall program
All the features of the obe and mbe routines are accessible through the driveall program, with the exceptions of two of the most specialised ones (setting collapse operators explicitly and varying the number of sub-steps inside each time step of a time-dependent integration).There are also minor restrictions on certain modes of operation, as flagged in the User Manual.
All the necessary data and control parameters are passed to driveall through several input files.A computation using this program simply involves (i) specifying the required number of states and other key parameters in the general settings module; (ii) compiling the program; (iii) preparing or updating the input files as necessary; and (iv) executing the program.The program can be compiled once and for all, as long as the number of states and other key parameters specified in the general settings module are kept the same.
The driveall program reads up to five different input files, of which two must always be provided.The two mandatory files are referred to as the keyparams file and the controlparams file.Their content is listed in tables 3, 4, 5 and 6.For convenience, parameters such as energies and dipole moments, which The Rabi frequencies Ω α;i j divided by 2π, in MHz, expressed as real numbers.
may be numerous in calculations on large multistate systems, can be specified in an auxiliary defaultdata file as an alternative to being listed in the controlparams file.The program reads these files using the namelist feature of Fortran, as explained in Appendix E and Appendix F. These files must therefore be formatted accordingly; however, no knowledge of Fortran programming is required beyond what is mentioned in this regard in that appendix.Further examples of the use of driveall can be found in the examples folder forming part of the distribution.The reader is referred to the User Manual for a full description of all the features of this program.

Running these codes through a bespoke program 4.5.1. Representation of the density matrix
As was explained in Section 3.1.1,the density matrices are stored within these modules as column vectors of N 2 real numbers, as per Eq. ( 17): a 1D array rhovec representing a density matrix is such that rhovec(1) contains ρ 11 , rhovec(2) contains Re ρ 12 , etc. (or ρ 00 , Re ρ 01 , etc., if the states are numbered 0, 1, 2,. . .rather than 1, 2, 3,. . ., see Section 4.6).Which components of such vectors correspond to which elements of the density matrix can be found by using the subroutines obe coher index and obe pop index of the obe module.

The obefield and obecfield derived types
Two Fortran derived variable types, called obefield and obecfield, are defined in the module obe.Variables of this type are used by the obe and mbe programs for storing and communicating various attributes of the relevant fields, such as their amplitude, wavelength, direction of propagation, detuning, and the Rabi frequencies or dipole moments of the transitions they drive.A full description of these two derived types can be found in the User Manual.No knowledge of these derived types is required for running the codes through the driveall program.

Structure of the program
With the exceptions of the subroutines flagged at the end of this section, using this package will normally involve the steps outlined below.A Fortran 90 program using the obe steadystate subroutine for a steady-state calculation is provided in Appendix E as an example.A program using the mbe module for a propagation calculation is also included in the examples directory.Detailed information about using the various user-facing routines provided in this library can be found in the User Manual.
1.The driving program should first pass various pieces of information to the obe module through a call to the subroutine obe setcsts, namely the frequency offsets of the different states, the rates of spontaneous decay and optionally any additional dephasing rate and any additional collapse operator, as well as the number of fields, whether the weak probe approximation is to be assumed or not, and whether Rabi frequencies or complex electric field amplitudes and dipole moments will be used for defining how each of the fields couples to the different states.
2. The parameters of each of the fields must then be passed to obe through calls to the subroutine obe setfields.
The field identified by the reference number 1 in the corresponding call to obe setfields is taken to be the probe field if the calculation is to be done within the weak probe approximation.
3. The root-mean squared velocity of the atoms and the details of the integration over atomic velocities must be passed to obe through a call to the subroutine obe set Doppler if a calculation involving a Doppler averaging by numerical quadrature is to be done.A choice of general numerical quadratures is offered by obe set Doppler.Alternatively, the user can upload custom abscissas and weights.
4. Unless the applied fields are CW, the details of their temporal envelope must be passed to mbe either through a call to mbe set envlp followed by a call to mbe set tdfields A, or through a call to mbe set tdfields B. The latter makes it possible to use time meshes and define temporal profiles more varied than offered by mbe set envlp and mbe set tdfields A.
5. The relative and absolute accuracy parameters of the DOP853 ODE solver must also be passed to obe, through a call to obe set tol dop853, if this solver is to be used in the course of the calculation.
6. Unit numbers for the output of selected elements of the density matrix must be passed to obe through a call to obe setoutputfiles if this option of outputting results is to be used.
7. The relevant computational routines must then be called for performing the required calculation.The initial populations and (possibly) coherences need to be passed to these various subroutines as input data, with the exceptions mentioned in their detailed descriptions in the User Manual.
8. A calculation of the density matrix can be followed, if required, by a call to obe susceptibility, which calculates the complex susceptibility, refractive index and absorption coefficient.9.For propagation calculations, the program must also include an external subroutine through which mbe propagate 1 or mbe propagate 2 can output the results, as described in the User Manual.
The detunings and complex amplitudes of the applied fields initially set by obe setfields can be reset at a later stage, respectively by calling the subroutines obe reset detuning and obe reset cfield of the obe module.This makes it possible, e.g., to calculate refractive indexes and absorption coefficients for a range of detunings or a range of field strengths.
The obe module also includes several auxiliary routines which may be of assistance when preparing the input of some of the subprograms mentioned above or processing their output.These are obe coher index and obe pop index, for identifying the relevant elements of a density matrix in the 1D storage mode described in Section 4.5.1;obe find cfield and obe find rabif, for relating complex electric field amplitudes to complex Rabi frequencies in the definition of Eq. ( 10); and obe init rho, for initialising a density matrix in its 1D representation.
None of the initialisation steps listed above are necessary if the only computational routines used would be obe 2state, obe weakprb 3stladder, obe weakprb 4stladder or obe weakfield.

State numbering
The N states included in a calculation are identified by numbers running from 1 to N throughout these code and in the documentation.This numbering is in line with the default indexing of arrays in Fortran.However, the user may choose to use a different numbering system for describing the system in the driving program, e.g., one where the states are identified by a number running from 0 to N −1 rather from 1 to N. The value of the variable nmn set in the general settings module informs the user-facing obe and mbe routines of the numbering system used in the external programs calling them: giving a value of n to nmn means that the states are numbered n, n + 1, n + 2,. . . in the information passed to obe and mbe by the driving program.E.g., setting nmn to 0 means that the states are numbered 0, 1, 2,. . . in the arrays passed to these modules by the user, while setting nmn to 1 means that the states are instead numbered 1, 2, 3,. . .Within the obe, mbe and ldbl modules, however, the states are numbered 1, 2 and 3, irrespective of the value of nmn.
For instance, the three states of the systems considered in Appendix E could be identified equally well by the numbers 0, 1 and 2, rather than by the numbers 1, 2 and 3.In order to use 0, 1 and 2, the constant nmn should be given the value 0 in the general settings module (and in the keyparams file if driveall is used).The statements Gamma_decay(1,2) = 5.0d0 Gamma_decay(2,3) = 1.0d0 defining the decay rates in the example given in Appendix E should then be replaced by Gamma_decay(0,1) = 5.0d0 Gamma_decay(1,2) = 1.0d0 and similarly for the arrays Rabif, detuning fact and energ f.

Code reuse
The ldbl and mbe modules both contain a copy, in essentially the original form, of the subroutine DOP853 described in Ref. [27] and published by the University of Geneva [28].The obe module contains a copy of the subroutine CLENSHAW CURTIS COMPUTE published by J. Burkardt [35], and a Fortran 90 implementation of the wwerf function of the CERN Library, which calculates the Faddeeva function [36,37,38].
The structure of the optical Bloch equations ensures that the L iJ 's forming the right-hand sides of Eq. (B.4) do not depend on detunings, and that the other L i j 's depend at most linearly on them.The column vector b is thus constant in the atomic velocity v, while the matrix L ′ varies linearly with v.We set, accordingly, where L ′ 0 and L ′ 1 do not depend on v.The two matrices L ′ 0 and L ′ 1 define the generalized eigenvalue problem where µ is a (normally complex) generalized eigenvalue.Since L ′ 0 and L ′ 1 are (N −1)×(N −1) matrices, the span of the solution vectors x is a space of dimension M ≤ N − 1.It is thus possible to find M eigenvectors x 1 , x 2 ,. . ., x M forming a basis for this space.With µ j denoting the corresponding eigenvalues, To each eigenvector x j can be associated a left eigenvector y j such that and The solution r ′ of Eq. ( 33) can be written as the sum of a linear combination of the eigenvectors x j 's and of a vector r ′ 0 biorthogonal to all the left eigenvectors: with r ′ 0 being such that (Formally, c j = y † j L ′ 1 r ′ and r ′ 0 = r ′ − j c j x j .)Combining the above equations yields and Solving Eqs.(B.8) and (B.9) for the eigenvalues µ j and the corresponding right and left eigenvectors is a standard numerical problem, as is solving Eq. (B.14) for the vector r ′ 0 .(In the present programs, this calculation is done by first reverting to a formulation of the density matrix in terms of real populations and complex coherences, and working with the complex matrices and complex vectors corresponding to L ′ , r ′ and b in that formulation.)Altogether, the calculation yields each of the elements of r st as a sum of partial fractions with constant numerators and denominators linear in v: where the α i j 's are constants.As is well known, expressions of this form are readily amenable to an analytical averaging over a Maxwellian distribution of velocities [29], and indeed, expanding coherences in partial fractions of this form is a standard approach in few-state calculations based on the weak probe approximation.The method outlined in this appendix generalises this approach to multi-state, multi-fields systems treated beyond the weak probe approximation.Doppler averaging is based on the following identities, where η j = −µ j /u and, as defined in Section 3.1.2,w(•) is the Faddeeva function: (B.17 The case Im η j = 0 does not need to be considered as the eigenvalues µ j always have a non-zero imaginary part for any pair of matrices L ′ 0 and L ′ 1 arising from the optical Bloch equations.Organising the calculations along similar lines may also lead to a significant speed up in computations of the steady state density matrix for multiple values of one of the detunings.For such calculations, Eq. (B.6) would be replaced by the equation where the matrices L′ 0 and L′ 1 do not depend on ∆ α .Following the above procedure then results in a density matrix of the form where αi j and μ j are constants.The only potentially CPU intensive step in this approach is the calculation of the generalized eigenvalues and eigenvectors of the matrix pair ( L′ 0 , L′ 1 ), which does not need to be repeated for each value of ∆ α .

Appendix C. The optical Bloch equations in the weak probe approximation
For simplicity, we will only consider the case where the amplitude of the probe field, E p , is real.The final results -Eqs.(C.6) and (C.7) -are easily generalised to the case of a complex amplitude, and the program is organised in such a way that the weak probe approximation is correctly implemented whether the amplitude of the probe field is real or complex.
We assume that the coherences are all initially zero.The populations then vary with E p only through terms quadratic or of higher order in E p .The populations will therefore vary little if the probe field is very weak.The essence of the weak probe approximation is to integrate Eq. ( 18) only to the leading (non vanishing) order in E p .This is done to first order in E p within the obe module.
Implementing this approximation first requires a consideration of Eq. ( 18) in the limit of a vanishing probe field (E p → 0).The elements of the density matrix divide into two classes in that limit, namely the populations and the coherences which take on non-zero values either initially or at later times (class A), and the coherences which are initially zero and remain zero at all times (class B).(The elements of class A may vary in time even when E p = 0, e.g., because of spontaneous decay or because of an interaction with a field other than the probe field.)We can thus form the column vector r by concatenating the column vectors formed by the respective populations and coherences, r A and r B : Eq. ( 15) simplifies considerably in that limit: ρ i j ≡ 0 if i and j both belong to G 1 or both belong to G 2 , whereas when i ∈ G 2 and j ∈ G 1 .The decoherence rates γ i j account for dephasing mechanisms not contributing to the decay rates Γ i , such as random phase jumps of the field contributing to its frequency width [39] and collisional broadening.Typically, where γ coll i j is the decay rate of the coherence ρ i j due to collisions and ∆ν is the frequency width of the field (full width at half maximum).
We refer the energies of the states to either an energy ℏω ref (1) or ℏω ref (2) depending on whether they belong to group G 1 or group G 2 .Thus and therefore, in the above equation, Moreover, since we defined G 2 as containing states higher in energy than the states belonging to G 1 , Setting ρi j = 0 yields the steady state coherences: with Given Eqs.(D.8) and (D.9), Eq. ( 36) yields a particularly simple result for the corresponding complex susceptibility: Note that the full width at half maximum of the resonance peak at ∆ i j = 0 is twice the total dephasing rate γ tot i j .E.g., to obtain a collisional width of Γ coll i j (full width at half maximum in angular frequency), the dephasing rate γ coll i j must be set equal to Γ coll i j /2.Doppler averaging χ(ω 1 ) then amounts to a simple application of Eqs.(B.16) and (B.17), since with η i j = ∆ i j + i γ tot i j /(uk).Therefore (D.13)

Appendix E. Example of steady state calculation
This appendix illustrates how the obe codes can be used for calculating the steady state density matrix for a ladder system of three states, states 1, 2 and 3, with ℏω (1) < ℏω (2) < ℏω (3) .States 1 and 2 are coupled to each other by field 1 (the "probe field") and states 2 and 3 by field 2 (the "coupling field").Within the rotating wave approximation, the Hamiltonian of this systems is represented by the following matrix in the {|1⟩, |2⟩, |3⟩} basis, where ∆ 1 and ∆ 2 are given by Eqs. ( 13) and (14).It is assumed that state 3 decays to state 2 and state 2 to state 1, the corresponding decay rates being Γ 32 and Γ 21 , respectively, and that inhomogeneous broadening can be neglected.Specifically, we take δω (1) = δω (2) = δω We first show how the steady state density matrix could be obtained by running the code through the driveall program.We then give an example of a bespoke program doing the same calculation.Copies of these files are provided in the examples folder.

Using the driveall program
Calculating this density matrix using the driveall program requires two input files, namely the keyparams and the controlparams files, formatted as illustrated by the examples below.Apart from possible comments and blank lines, each of these two files must start with an ampersand symbol followed by the name of the respective namelist structure (keyparams for the keyparams file, controlparams for the controlparams file) and must end with a slash.Each input value must be provided in the form of a Fortran assignment statement (e.g., i = 1 if i is an integer variable, v = 1.0d0 if v is a double precision variable, s = 'something' if s is a character variable).Input values can be provided in any order and do not need to be all present.Strings of characters starting with an exclamation mark are taken to be comments and are ignored, as are blank lines.
The following file could be used as the keyparams file for that calculation.This file is read by driveall from the standard input stream.It specifies several key parameters and the name of the controlparams file.!Number of fields.icmplxfld = 0 !Indicates that the field amplitudes ! and Rabi frequencies will be specified as real !numbers, not as complex numbers.
! Name of the controlparams files: filename_controlparams = 'example_c.dat'/ The corresponding controlparams file could be taken to be as follows. 3controlparams icalc = 2 !Tells driveall to calculate the !steady state.iRabi = 1 !The Rabi frequencies will be specified.inoncw = 0 !The fields are CW.iweakprb = 0 !The weak probe approximation is not made.iDoppler = 0 !No Doppler averaging.!Rabi frequencies, in units of ( 2 !The detunings, in units of (2pi) x MHz: detuning(1) = 5.0d0 detuning(2) = 0.0d0 / Running driveall with these input files produces the following output:   [42].The spontaneous decay rates of the 5 2 P 1/2 and 5 2 P 3/2 states are, respectively, 2π × 5.746 MHz and 2π × 6.0666 MHz for the 5 2 P 3/2 state [16,41].The atomic number density, 1.96 × 10 21 m −3 , corresponds to a vapour temperature of 220 • C; however, Doppler broadening is neglected in this example so as to avoid unnecessarily long execution times.The Maxwell-Bloch equations are integrated in time using Butcher's 5th order Runge-Kutta method and in space using a predictor-corrector method combining the third order Adams-Bashford method and the fourth order Adams-Moulton method, initiated by calculation with smaller spatial steps using the classic fourth order Runge-Kutta rule.The applied fields are read from a file called appliedfields.dat.This file, as well as the controlparams file listed below and the corresponding keyparams file can be found in the examples directory included in this distribution.Running the driveall program with these input data produces a file called outamplitudes.dat containing the complex amplitudes of the propagated fields as functions of position and time (however, see Appendix G for adapting the program listed below to being run through Podman).These complex amplitudes are transformed into the corresponding intensities by the program used for plotting Fig.The examples directory also contains a bespoke program doing the same calculation, although with the fields calculated directly within the mbe module rather than read from file.

Appendix G. Running CoOMBE using a container image
This code can be compiled and run using a container image such as one managed by the Podman tool [23].An advantage of this method is that the user does not need to worry about installing a Fortran compiler or the required numerical libraries.This approach is becoming increasingly conventional across modern software development.The user wishing to use Podman will need to install this program for their respective operating system (this program is freely available at [23]).Once done, the user can then easily build the relevant image of the code and run the program.
We provide a Podman implementation for each of the worked examples in the examples folder, namely a Containerfile (here a Dockerfile), a Makefile and a .shshell script.Building the image is done by the following command: podman build -t coombe .
(The ldbl.f90,obe.f90, mbe.f90 and driveall.f90files need to be first copied into the working directory as necessary, together with the relevant general settings.f90file, data files, Dockerfile, Makefile and shell script.)Once the image has been built, the program can be run using the command podman run -v ./:/home/coombecoombe Changing any input parameters normally requires to rebuild the image by using the podman build command again; however, the rebuild process is typically faster than the initial rebuild.
Table 5: Contents of the controlparams file read by the driveall program: II.Parameters specific to particular types of calculation.Comprehensive information about these parameters and the use of driveall can be found in the User Manual.The parameters indicated by an asterisk may be specified in the defaultdata file rather than in the controlparams file.

Name
Short description Parameters to be provided for time-dependent calculations of the density matrix imethod Parameter determining which numerical algorithm is to be used in the integration of the optical Bloch equation.n time steps N t , the number of integration steps in the time integration of the optical Bloch equations.(*) popinit The initial populations.
Optional parameters relevant for time-dependent calculations of the density matrix filename rhoall out The name of the file to which the whole density should be written at each time step.iAorB Parameter determining which subroutine should be used in calculations with Doppler averaging.inoncw Parameter determining whether the calculation is for CW fields or for fields with a time-dependent envelope.iprintrho Parameter determining whether the final density matrix is to be written out to the standard output stream.irate Parameter determining whether the calculation is to be done within the rate equations approximation.iunformatted Parameter determining whether the output files should be unformatted (binary) rather than formatted.

ti, tf
The lower and upper bounds of the time integration interval, in µs.
Optional parameters relevant for calculations of the steady state density matrix filename chi out The name of the file to which the calculated complex susceptibilities, refractive indexes and absorption coefficients should be written.iDoppler numer st Parameter determining whether Doppler averaging is to be done analytically or numerically.iladder wkprb Parameter determining whether the steady state density matrix is to be calculated by a subroutine specialised to ladder systems rather than by general subroutines.ioption Parameter determining the algorithm used for calculating the steady state density matrix.iprintrho Parameter determining whether the steady state density matrix is to be written out to the standard output stream.isuscept Parameter determining whether the complex susceptibility at the probe frequency and the corresponding values of the refractive index and absorption coefficient are calculated after the steady state density matrix has been obtained.ivarydetuning Parameter determining whether the steady state density matrix is to be calculated over a range of detunings.(*) popinit The initial populations.
Parameters to be provided in the case of propagation calculations (*) density The density of the medium expressed as the number of atoms per m 3 .imethod Parameter determining which numerical method is to be used for integrating the optical Bloch equations.n time steps N t , the number of integration steps in the time integration of the optical Bloch equations.n z steps The number of integration steps to be taken between z = 0 and z = z max .(*) popinit The initial populations.wavelength The wavelength(s) of the field(s) considered in nm.zmax The distance over which the field(s) must be propagated, z max , in µm.

Optional parameters relevant for propagation calculations filename rhoall out
The name of the file to which the whole density matrix should be written.iunformatted Parameter determining whether the output files should be unformatted (binary) rather than formatted.izrule Parameter determining the numerical algorithm used in the spatial propagation.nt writeout, nz writeout Constants determining at which time or spatial steps results should be written out.

ti, tf
The lower and upper bounds of the time integration interval, in µs.
Table 6: Contents of the controlparams file read by the driveall program: III.Miscellaneous parameters.Comprehensive information about these parameters and the use of driveall can be found in the User Manual.The parameters indicated by an asterisk may be specified in the defaultdata file rather than in the controlparams file.

Name Short description
Parameters to be provided for calculations with analytical Doppler-averaging idir Parameter(s) indicating whether the corresponding fields propagate in the positive z-direction or the negative zdirection.(*) urms The root-mean squared velocity of the atoms in the laser propagation direction, u, in m s −1 .(*) wavelength The wavelength(s) of the field(s) considered, in nm.
Parameters to be provided for calculations with numerical Doppler-averaging filename Dopplerquad The name of the file containing the quadrature abscissas and quadrature weights to be used in the calculation, if this file is required.idir Parameter(s) indicating whether the corresponding fields propagate in the positive z-direction or the negative zdirection.irule Parameter determining which quadrature rule is to be used in the calculation.

n v values
The number of integration points in the integration over v.

(*) urms
The root-mean squared velocity of the atoms in the laser propagation direction, u, in m s −1 .vmax The maximum value of |v| to be considered, in m s −1 .(*) wavelength The wavelength(s) of the field(s) considered, in nm.
Parameters to be provided when the DOP853 ODE integrator is to be used The wavelength of the probe field in nm.

F. 1 ,
which is also included in the examples directory.

atol
Parameter controlling the allowed absolute error on the populations and coherences calculated by the program.rtol Parameter controlling the allowed relative error on the populations and coherences calculated by the program.Parameters relevant for calculations involving non-CW fields filename tdamps in The name of the file containing the time-dependent amplitude(s) of the field(s) considered, if this file is required.filename tdamps out The name of the file to which the program should write the time-dependent amplitude(s) of the field(s) considered.iforce0 Parameter(s) determining whether the corresponding field should be taken to be initially zero.iinterp Parameter determining whether tabulated field amplitudes should be interpolated.istart Parameter determining the initial values of the populations and coherences.itdfieldsAorB Parameter determining whether the amplitude(s) of the applied field(s) must be calculated by the program or read from file.nsubsteps Number of substeps within each integration step.pulse type Parameter(s) determining the pulse envelope of the respective field.t0, t1, tw Parameters determining the shape of the pulse envelope of the respective field.Parameters to be provided for steady state calculations over a range of detunings detuning min, detuning max The smallest and largest values of the detuning divided by 2π, in MHz.detuning step The step in detuning, divided by 2π, in MHz.index field Parameter determining which field should be varied in the calculation.Parameters to be provided for a calculation of the susceptibility at the probe frequency (*) density The density of the medium expressed as the number of atoms per m 3 .(*) wavelength

Table 1 :
Contents of the distribution

Table 2 :
User-facing subprograms contained in the obe and mbe modules.For a single field with a time-dependent envelope, with or without Doppler averaging.

Table 3 :
Contents of the keyparams file read by the driveall program.Comprehensive information about these parameters and the use of driveall can be found in the User Manual.

Table 4 :
Contents of the controlparams file read by the driveall program: I. General parameters.Comprehensive information about these parameters and the use of driveall can be found in the User Manual.The parameters indicated by an asterisk may be specified in the defaultdata file rather than in the controlparams file.The electric field amplitudes E α in V m −1 , expressed as real numbers.(*)camplitudeTheelectricfieldamplitudes E α in V m −1 , expressed as complex numbers.(*)cdipmomTheelectricdipole moments⟨ i | εα • D | j ⟩ in C m,expressed as complex numbers.(*)cRabif The Rabi frequencies Ω α;i j divided by 2π, in MHz, expressed as complex numbers.(*)detuning The detunings ∆ α divided by 2π, in MHz.(*) detuning fact The detuning factors a iα .(*)dip mom The electric dipole moments ⟨ i | εα • D | j ⟩ in C m, expressed as real numbers.(*)energ f The energy offsets δω (i) divided by ℏ, in MHz.filename rho out The name(s) of the output file(s) to which the program should write the density matrix.(*)Gamma decay f The spontaneous decay rates Γ i j divided by 2π, in MHz.iappend Parameter determining whether existing output files can be overwritten with new results.iDopplerParameter determining whether the density matrix must be Doppler-averaged.iweakprbParameterdetermining whether the calculation is to be done within the weak probe approximation.
(*) Rabif The program listed below calculates and writes out the steady state value of ρ 12 for the same parameters.Running it produces the following output: