Notebook[{ Cell[CellGroupData[{ Cell["Find evolution equations", "Section"], Cell["Load the package.", "Text"], Cell[BoxData[ RowBox[{"<<", "AtomicDensityMatrix`"}]], "Input", CellGroupingRules->{GroupTogetherGrouping, 10001.}, CellID->2058623809], Cell["Use density matrix variables with explict time dependence.", "Text"], Cell[BoxData[ RowBox[{"SetOptions", "[", RowBox[{"DensityMatrix", ",", RowBox[{"TimeDependence", "\[Rule]", "True"}]}], "]"}]], "Input"], Cell["\<\ Pull quantum numbers and other basic info about the states from a database.\ \>", "Text", CellID->429217524], Cell[BoxData[ RowBox[{"s12data", "=", RowBox[{"AtomicData", "[", RowBox[{ RowBox[{"{", RowBox[{"\"\\"", ",", "87", ",", RowBox[{"{", RowBox[{"\"\\"", ",", RowBox[{"{", RowBox[{"5", ",", "\"\\""}], "}"}]}], "}"}], ",", RowBox[{"{", RowBox[{"2", ",", "\"\\"", ",", FractionBox["1", "2"]}], "}"}]}], "}"}], ",", RowBox[{"{", RowBox[{ "Energy", ",", "J", ",", "L", ",", "S", ",", "NuclearSpin", ",", "NaturalWidth", ",", "Parity"}], "}"}]}], "]"}]}]], "Input"], Cell[BoxData[ RowBox[{"p12data", "=", RowBox[{"Append", "[", RowBox[{ RowBox[{"AtomicData", "[", RowBox[{ RowBox[{"{", RowBox[{"\"\\"", ",", "87", ",", RowBox[{"{", RowBox[{"\"\\"", ",", RowBox[{"{", RowBox[{"5", ",", "\"\\""}], "}"}]}], "}"}], ",", RowBox[{"{", RowBox[{"2", ",", "\"\\"", ",", FractionBox["1", "2"]}], "}"}]}], "}"}], ",", RowBox[{"{", RowBox[{"J", ",", "L", ",", "S", ",", "NuclearSpin", ",", "Parity"}], "}"}]}], "]"}], ",", RowBox[{ RowBox[{"BranchingRatio", "[", "0", "]"}], "\[Rule]", "1"}]}], "]"}]}]], "Input"], Cell[BoxData[ RowBox[{"p32data", "=", RowBox[{"Append", "[", RowBox[{ RowBox[{"AtomicData", "[", RowBox[{ RowBox[{"{", RowBox[{"\"\\"", ",", "87", ",", RowBox[{"{", RowBox[{"\"\\"", ",", RowBox[{"{", RowBox[{"5", ",", "\"\\""}], "}"}]}], "}"}], ",", RowBox[{"{", RowBox[{"2", ",", "\"\\"", ",", FractionBox["3", "2"]}], "}"}]}], "}"}], ",", RowBox[{"{", RowBox[{"J", ",", "L", ",", "S", ",", "NuclearSpin", ",", "Parity"}], "}"}]}], "]"}], ",", RowBox[{ RowBox[{"BranchingRatio", "[", "0", "]"}], "\[Rule]", "1"}]}], "]"}]}]], "Input"], Cell["\<\ Using the quantum numbers, create a list of all hyperfine and Zeeman \ sublevels of the J states, which are labeled 0, 1, and 2 for reference. This \ list will be passed to the functions from the ADM package that create the DM \ evolution equations.\ \>", "Text"], Cell[BoxData[ RowBox[{ RowBox[{"system0", "=", RowBox[{"Sublevels", "[", RowBox[{"{", RowBox[{ RowBox[{"AtomicState", "[", RowBox[{"0", ",", "s12data"}], "]"}], ",", RowBox[{"AtomicState", "[", RowBox[{"1", ",", "p12data"}], "]"}], ",", RowBox[{"AtomicState", "[", RowBox[{"2", ",", "p32data"}], "]"}]}], "}"}], "]"}]}], ";"}]], "Input", CellID->433132487], Cell[TextData[{ "For simplicity, we delete some excited states from the system, leaving only \ ", Cell[BoxData[ FormBox[ RowBox[{ RowBox[{ SubscriptBox["P", RowBox[{"1", "/", "2"}]], " ", "F"}], "=", "1"}], TraditionalForm]]], " and ", Cell[BoxData[ FormBox[ RowBox[{ RowBox[{ SubscriptBox["P", RowBox[{"3", "/", "2"}]], " ", "F"}], "=", "2"}], TraditionalForm]]], "." }], "Text"], Cell[BoxData[ RowBox[{"system", "=", RowBox[{"DeleteStates", "[", RowBox[{"system0", ",", RowBox[{ RowBox[{"(", RowBox[{ RowBox[{"Label", "\[Equal]", "1"}], "&&", RowBox[{"F", "\[NotEqual]", "1"}]}], ")"}], "||", RowBox[{"(", RowBox[{ RowBox[{"Label", "\[Equal]", "2"}], "&&", RowBox[{"F", "\[NotEqual]", "2"}]}], ")"}]}]}], "]"}]}]], "Input"], Cell[TextData[{ "Define the optical field with four frequencies, ", Cell[BoxData[ FormBox[ StyleBox[ SubscriptBox["\[Omega]", "1"], "InlineMath"], TraditionalForm]]], ", ", Cell[BoxData[ FormBox[ StyleBox[ SubscriptBox["\[Omega]", "2"], "InlineMath"], TraditionalForm]]], ", ", Cell[BoxData[ FormBox[ StyleBox[ SubscriptBox["\[Omega]", "3"], "InlineMath"], TraditionalForm]]], ", and ", Cell[BoxData[ FormBox[ StyleBox[ SubscriptBox["\[Omega]", "4"], "InlineMath"], TraditionalForm]]], ". Here ", Cell[BoxData[ FormBox[ SubscriptBox["\[CapitalOmega]", "1"], TraditionalForm]]], " and ", Cell[BoxData[ FormBox[ SubscriptBox["\[CapitalOmega]", "1"], TraditionalForm]]], " are left-circularly polarized, and ", Cell[BoxData[ FormBox[ SubscriptBox["\[CapitalOmega]", "3"], TraditionalForm]]], " and ", Cell[BoxData[ FormBox[ SubscriptBox["\[CapitalOmega]", "4"], TraditionalForm]]], " are right-circularly polarized." }], "Text", CellID->133602844], Cell[BoxData[ RowBox[{"SetOptions", "[", RowBox[{"OpticalField", ",", RowBox[{"PolarizationVector", "\[Rule]", RowBox[{"{", RowBox[{"1", ",", "0", ",", "0"}], "}"}]}], ",", RowBox[{"CartesianCoordinates", "\[Rule]", RowBox[{"{", RowBox[{"x", ",", "y", ",", "z"}], "}"}]}]}], "]"}]], "Input"], Cell[TextData[{ "Parameters for the OpticalField function are {frequency, wavenumber}, \ {electric amplitude, phase}, {rotation angle (relative to \ PolarizationVector), ellipticity}. Here ", Cell[BoxData[ FormBox[ SubscriptBox["\[CapitalOmega]", "i"], TraditionalForm]]], " are the Rabi frequencies defined in terms of the dipole reduced matrix \ elements." }], "Text"], Cell[BoxData[ RowBox[{"field", "=", RowBox[{ RowBox[{"OpticalField", "[", RowBox[{ RowBox[{"{", RowBox[{ SubscriptBox["\[Omega]", "1"], ",", SubscriptBox["k", "1"]}], "}"}], ",", RowBox[{"{", RowBox[{ RowBox[{ SubscriptBox["\[CapitalOmega]", "1"], "/", RowBox[{"ReducedME", "[", RowBox[{"0", ",", RowBox[{"{", RowBox[{"Dipole", ",", "1"}], "}"}], ",", "1"}], "]"}]}], ",", SubscriptBox["\[Phi]", "1"]}], "}"}], ",", RowBox[{"{", RowBox[{"0", ",", RowBox[{"\[Pi]", "/", "4"}]}], "}"}]}], "]"}], "+", "\[IndentingNewLine]", RowBox[{"OpticalField", "[", RowBox[{ RowBox[{"{", RowBox[{ SubscriptBox["\[Omega]", "2"], ",", SubscriptBox["k", "2"]}], "}"}], ",", RowBox[{"{", RowBox[{ RowBox[{ SubscriptBox["\[CapitalOmega]", "2"], "/", RowBox[{"ReducedME", "[", RowBox[{"0", ",", RowBox[{"{", RowBox[{"Dipole", ",", "1"}], "}"}], ",", "1"}], "]"}]}], ",", SubscriptBox["\[Phi]", "2"]}], "}"}], ",", RowBox[{"{", RowBox[{"0", ",", RowBox[{"\[Pi]", "/", "4"}]}], "}"}]}], "]"}], "+", "\[IndentingNewLine]", RowBox[{"OpticalField", "[", RowBox[{ RowBox[{"{", RowBox[{ SubscriptBox["\[Omega]", "3"], ",", SubscriptBox["k", "3"]}], "}"}], ",", RowBox[{"{", RowBox[{ RowBox[{ SubscriptBox["\[CapitalOmega]", "3"], "/", RowBox[{"ReducedME", "[", RowBox[{"0", ",", RowBox[{"{", RowBox[{"Dipole", ",", "1"}], "}"}], ",", "2"}], "]"}]}], ",", SubscriptBox["\[Phi]", "3"]}], "}"}], ",", RowBox[{"{", RowBox[{"0", ",", RowBox[{ RowBox[{"-", "\[Pi]"}], "/", "4"}]}], "}"}]}], "]"}], "+", "\[IndentingNewLine]", RowBox[{"OpticalField", "[", RowBox[{ RowBox[{"{", RowBox[{ SubscriptBox["\[Omega]", "4"], ",", SubscriptBox["k", "4"]}], "}"}], ",", RowBox[{"{", RowBox[{ RowBox[{ SubscriptBox["\[CapitalOmega]", "4"], "/", RowBox[{"ReducedME", "[", RowBox[{"0", ",", RowBox[{"{", RowBox[{"Dipole", ",", "1"}], "}"}], ",", "2"}], "]"}]}], ",", SubscriptBox["\[Phi]", "4"]}], "}"}], ",", RowBox[{"{", RowBox[{"0", ",", RowBox[{ RowBox[{"-", "\[Pi]"}], "/", "4"}]}], "}"}]}], "]"}]}]}]], "Input", CellID->534530029], Cell["\<\ The Hamiltonian for the system subject to the optical field. Each field is \ assumed to interact with only one transition\[LongDash]the replacement rule \ (Cos[_]|Sin[_]) ReducedME[_,{Dipole,1},_]\[Rule]0 causes other terms to be \ set to zero.\ \>", "Text", CellID->462076121], Cell[BoxData[ RowBox[{"MatrixForm", "[", RowBox[{"H", "=", RowBox[{ RowBox[{"Expand", "@", RowBox[{"Hamiltonian", "[", RowBox[{"system", ",", RowBox[{"ElectricField", "\[Rule]", "field"}], ",", RowBox[{"MagneticField", "\[Rule]", RowBox[{"{", RowBox[{"0", ",", "0", ",", RowBox[{"\[CapitalOmega]L", "/", "BohrMagneton"}]}], "}"}]}]}], "]"}]}], "/.", " ", RowBox[{ RowBox[{ RowBox[{"(", RowBox[{ RowBox[{"Cos", "[", "_", "]"}], "|", RowBox[{"Sin", "[", "_", "]"}]}], ")"}], " ", RowBox[{"ReducedME", "[", RowBox[{"_", ",", RowBox[{"{", RowBox[{"Dipole", ",", "1"}], "}"}], ",", "_"}], "]"}]}], "\[Rule]", "0"}]}]}], "]"}]], "Input", CellID->494599775], Cell["\<\ The level diagram for the system, showing optical couplings. Note that both \ resonant and off-resonant (counter-rotating) couplings are shown, because we \ have not yet performed the rotating-wave approximation.\ \>", "Text", CellID->358620443], Cell[BoxData[ RowBox[{"LevelDiagram", "[", RowBox[{"system", ",", RowBox[{"H", "/.", RowBox[{"{", RowBox[{ RowBox[{ RowBox[{"Energy", "[", "0", "]"}], "\[Rule]", "0"}], ",", RowBox[{ RowBox[{"Energy", "[", "1", "]"}], "\[Rule]", "2"}], ",", RowBox[{ RowBox[{"Energy", "[", "2", "]"}], "\[Rule]", "4"}], ",", RowBox[{ RowBox[{"HyperfineA", "[", "_", "]"}], "\[Rule]", ".2"}], ",", RowBox[{ RowBox[{"HyperfineB", "[", "_", "]"}], "\[Rule]", ".1"}], ",", RowBox[{"\[CapitalOmega]L", "\[Rule]", ".1"}]}], "}"}]}], ",", RowBox[{"ParityOffset", "\[Rule]", "False"}]}], "]"}]], "Input", CellID->167259034], Cell[TextData[{ "Here we apply the rotating-wave approximation to the Hamiltonian. Here we \ construct a list of frequency shifts for the list of Zeeman sublevels that \ will have the effect of removing the optical frequencies from the \ Hamiltonian. I.e., we hold ", Cell[BoxData[ FormBox[ RowBox[{ RowBox[{ SubscriptBox["S", RowBox[{"1", "/", "2"}]], "F"}], "=", "2"}], TraditionalForm]]], " fixed, shift ", Cell[BoxData[ FormBox[ SubscriptBox["P", RowBox[{"1", "/", "2"}]], TraditionalForm]]], " down by ", Cell[BoxData[ FormBox[ SubscriptBox["\[Omega]", "1"], TraditionalForm]]], ", shift ", Cell[BoxData[ FormBox[ RowBox[{ RowBox[{ SubscriptBox["S", RowBox[{"1", "/", "2"}]], " ", "F"}], "=", "1"}], TraditionalForm]]], " down by ", Cell[BoxData[ FormBox[ RowBox[{ SubscriptBox["\[Omega]", "1"], "-", SubscriptBox["\[Omega]", "2"]}], TraditionalForm]]], ", and shift ", Cell[BoxData[ FormBox[ SubscriptBox["P", RowBox[{"3", "/", "2"}]], TraditionalForm]]], " down by ", Cell[BoxData[ FormBox[ RowBox[{ SubscriptBox["\[Omega]", "1"], "-", SubscriptBox["\[Omega]", "2"], "+", SubscriptBox["\[Omega]", "3"]}], TraditionalForm]]], "." }], "Text", CellID->577766068], Cell[BoxData[{ RowBox[{ RowBox[{ RowBox[{"{", RowBox[{ RowBox[{"Label", "[", "#", "]"}], ",", RowBox[{"F", "[", "#", "]"}]}], "}"}], "&"}], "/@", "system"}], "\[IndentingNewLine]", RowBox[{"shifts", "=", RowBox[{"%", "/.", RowBox[{"{", RowBox[{ RowBox[{ RowBox[{"{", RowBox[{"0", ",", "2"}], "}"}], "\[Rule]", "0"}], ",", RowBox[{ RowBox[{"{", RowBox[{"1", ",", "_"}], "}"}], "\[Rule]", SubscriptBox["\[Omega]", "1"]}], ",", RowBox[{ RowBox[{"{", RowBox[{"0", ",", "1"}], "}"}], "\[Rule]", RowBox[{ SubscriptBox["\[Omega]", "1"], "-", SubscriptBox["\[Omega]", "2"]}]}], ",", RowBox[{ RowBox[{"{", RowBox[{"2", ",", "_"}], "}"}], "\[Rule]", RowBox[{ SubscriptBox["\[Omega]", "1"], "-", SubscriptBox["\[Omega]", "2"], "+", SubscriptBox["\[Omega]", "3"]}]}]}], "}"}]}]}]}], "Input"], Cell["\<\ Make a unitary transformation matrix to implement the chosen frequency shifts.\ \>", "Text"], Cell[BoxData[ RowBox[{"MatrixForm", "[", RowBox[{"transmat", "=", RowBox[{"MatrixExp", "[", RowBox[{ RowBox[{ RowBox[{"-", "\[ImaginaryI]"}], " ", "t", " ", RowBox[{"DiagonalMatrix", "[", "shifts", "]"}]}], "+", RowBox[{"\[ImaginaryI]", " ", "z", " ", RowBox[{"DiagonalMatrix", "[", RowBox[{"shifts", "/.", RowBox[{"\[Omega]", "\[Rule]", "k"}]}], "]"}]}]}], "]"}]}], "]"}]], "Input"], Cell[TextData[{ "Write ", Cell[BoxData[ FormBox[ SubscriptBox["\[Omega]", "4"], TraditionalForm]]], " and ", Cell[BoxData[ FormBox[ SubscriptBox["k", "4"], TraditionalForm]]], " in terms of the non-degenerate detuning and phase-mismatch parameters \ \[Delta]\[Omega] and \[Delta]k, which we set here to zero for simplicity. We \ then apply the transform matrix to the Hamiltonian and set off resonant terms \ oscillating at harmonics of ", Cell[BoxData[ FormBox[ SubscriptBox["\[Omega]", "1"], TraditionalForm]]], ", ", Cell[BoxData[ FormBox[ SubscriptBox["\[Omega]", "2"], TraditionalForm]]], ", and ", Cell[BoxData[ FormBox[ SubscriptBox["\[Omega]", "3"], TraditionalForm]]], " to zero." }], "Text"], Cell[BoxData[ RowBox[{"MatrixForm", "[", RowBox[{"H1", "=", RowBox[{"RotatingWaveApproximation", "[", RowBox[{"system", ",", RowBox[{ RowBox[{"H", "/.", RowBox[{"{", RowBox[{ RowBox[{ SubscriptBox["\[Omega]", "4"], "\[Rule]", RowBox[{ SubscriptBox["\[Omega]", "1"], "-", SubscriptBox["\[Omega]", "2"], "+", SubscriptBox["\[Omega]", "3"], "-", "\[Delta]\[Omega]"}]}], ",", RowBox[{ SubscriptBox["k", "4"], "\[Rule]", RowBox[{ SubscriptBox["k", "1"], "-", SubscriptBox["k", "2"], "+", SubscriptBox["k", "3"], "-", "\[Delta]k"}]}]}], "}"}]}], "/.", RowBox[{"{", RowBox[{ RowBox[{"\[Delta]\[Omega]", "\[Rule]", "0"}], ",", RowBox[{"\[Delta]k", "\[Rule]", "0"}]}], "}"}]}], ",", RowBox[{"{", RowBox[{ SubscriptBox["\[Omega]", "1"], ",", SubscriptBox["\[Omega]", "2"], ",", SubscriptBox["\[Omega]", "3"]}], "}"}], ",", RowBox[{"TransformMatrix", "\[Rule]", "transmat"}]}], "]"}]}], "]"}]], "Input"], Cell[TextData[{ "There are remaining fast-oscillating terms in the Hamiltonian at the \ difference frequency between ", Cell[BoxData[ FormBox[ SubscriptBox["\[Omega]", "1"], TraditionalForm]]], " and ", Cell[BoxData[ FormBox[ SubscriptBox["\[Omega]", "2"], TraditionalForm]]], ", due to the fact that each optical field can interact with both \ ground-state hyperfine levels. We set these terms to zero, which is \ equivalent to assuming that each field interacts with only the ", Cell[BoxData[ FormBox[ RowBox[{"F", "=", "1"}], TraditionalForm]]], " or the ", Cell[BoxData[ FormBox[ RowBox[{"F", "=", "2"}], TraditionalForm]]], " ground-state hyperfine sublevel." }], "Text"], Cell[BoxData[ RowBox[{"MatrixForm", "[", RowBox[{"H2", "=", RowBox[{"H1", "/.", RowBox[{"{", RowBox[{ RowBox[{ SuperscriptBox["\[ExponentialE]", RowBox[{"\[Phi]_", "-", RowBox[{"\[ImaginaryI]", " ", "t", " ", SubscriptBox["\[Omega]", "1"]}], "+", RowBox[{"\[ImaginaryI]", " ", "t", " ", SubscriptBox["\[Omega]", "2"]}]}]], "\[Rule]", "0"}], ",", RowBox[{ SuperscriptBox["\[ExponentialE]", RowBox[{"\[Phi]_", "+", RowBox[{"\[ImaginaryI]", " ", "t", " ", SubscriptBox["\[Omega]", "1"]}], "-", RowBox[{"\[ImaginaryI]", " ", "t", " ", SubscriptBox["\[Omega]", "2"]}]}]], "\[Rule]", "0"}]}], "}"}]}]}], "]"}]], "Input"], Cell["\<\ Find the frequencies at which each field is assumed to be resonant with its F\ \[Rule]F\[CloseCurlyQuote] transition.\ \>", "Text"], Cell[BoxData[ RowBox[{"\[Omega]1res", "=", RowBox[{ RowBox[{ RowBox[{"Hamiltonian", "[", RowBox[{"{", RowBox[{"SelectState", "[", RowBox[{"system", ",", " ", RowBox[{ RowBox[{"Label", "\[Equal]", "1"}], "&&", RowBox[{"F", "\[Equal]", "1"}]}]}], "]"}], "}"}], "]"}], "\[LeftDoubleBracket]", RowBox[{"1", ",", "1"}], "\[RightDoubleBracket]"}], "-", RowBox[{ RowBox[{"Hamiltonian", "[", RowBox[{"{", RowBox[{"SelectState", "[", RowBox[{"system", ",", " ", RowBox[{ RowBox[{"Label", "\[Equal]", "0"}], "&&", RowBox[{"F", "\[Equal]", "2"}]}]}], "]"}], "}"}], "]"}], "\[LeftDoubleBracket]", RowBox[{"1", ",", "1"}], "\[RightDoubleBracket]"}]}]}]], "Input"], Cell[BoxData[ RowBox[{"\[Omega]2res", "=", RowBox[{ RowBox[{ RowBox[{"Hamiltonian", "[", RowBox[{"{", RowBox[{"SelectState", "[", RowBox[{"system", ",", " ", RowBox[{ RowBox[{"Label", "\[Equal]", "1"}], "&&", RowBox[{"F", "\[Equal]", "1"}]}]}], "]"}], "}"}], "]"}], "\[LeftDoubleBracket]", RowBox[{"1", ",", "1"}], "\[RightDoubleBracket]"}], "-", RowBox[{ RowBox[{"Hamiltonian", "[", RowBox[{"{", RowBox[{"SelectState", "[", RowBox[{"system", ",", " ", RowBox[{ RowBox[{"Label", "\[Equal]", "0"}], "&&", RowBox[{"F", "\[Equal]", "1"}]}]}], "]"}], "}"}], "]"}], "\[LeftDoubleBracket]", RowBox[{"1", ",", "1"}], "\[RightDoubleBracket]"}]}]}]], "Input"], Cell[BoxData[ RowBox[{"\[Omega]3res", "=", RowBox[{ RowBox[{ RowBox[{"Hamiltonian", "[", RowBox[{"{", RowBox[{"SelectState", "[", RowBox[{"system", ",", " ", RowBox[{ RowBox[{"Label", "\[Equal]", "2"}], "&&", RowBox[{"F", "\[Equal]", "2"}]}]}], "]"}], "}"}], "]"}], "\[LeftDoubleBracket]", RowBox[{"1", ",", "1"}], "\[RightDoubleBracket]"}], "-", RowBox[{ RowBox[{"Hamiltonian", "[", RowBox[{"{", RowBox[{"SelectState", "[", RowBox[{"system", ",", " ", RowBox[{ RowBox[{"Label", "\[Equal]", "0"}], "&&", RowBox[{"F", "\[Equal]", "1"}]}]}], "]"}], "}"}], "]"}], "\[LeftDoubleBracket]", RowBox[{"1", ",", "1"}], "\[RightDoubleBracket]"}]}]}]], "Input"], Cell[TextData[{ "Rewrite the frequencies ", Cell[BoxData[ FormBox[ SubscriptBox["\[Omega]", "1"], TraditionalForm]]], ", ", Cell[BoxData[ FormBox[ SubscriptBox["\[Omega]", "2"], TraditionalForm]]], ", and ", Cell[BoxData[ FormBox[ SubscriptBox["\[Omega]", "3"], TraditionalForm]]], " in terms of a detuning ", Cell[BoxData[ FormBox[ SubscriptBox["\[Delta]", "1"], TraditionalForm]]], ", ", Cell[BoxData[ FormBox[ SubscriptBox["\[Delta]", "2"], TraditionalForm]]], ", or ", Cell[BoxData[ FormBox[ SubscriptBox["\[Delta]", "3"], TraditionalForm]]], " from resonance with the appropriate transition between ground-state and \ excited state hyperfine levels, and subtract a constant term off of the \ diagonal to simplify the appearance." }], "Text"], Cell[BoxData[ RowBox[{"MatrixForm", "[", RowBox[{"Hrwa", "=", RowBox[{ RowBox[{"(", RowBox[{"H2", "/.", RowBox[{"{", RowBox[{ RowBox[{ SubscriptBox["\[Omega]", "1"], "\[Rule]", RowBox[{"\[Omega]1res", "+", SubscriptBox["\[Delta]", "1"]}]}], ",", RowBox[{ SubscriptBox["\[Omega]", "2"], "\[Rule]", RowBox[{"\[Omega]2res", "+", SubscriptBox["\[Delta]", "2"]}]}], ",", RowBox[{ SubscriptBox["\[Omega]", "3"], "\[Rule]", RowBox[{"\[Omega]3res", "+", SubscriptBox["\[Delta]", "3"]}]}]}], "}"}]}], ")"}], "-", RowBox[{ FractionBox[ RowBox[{"3", " ", RowBox[{"HyperfineA", "[", "0", "]"}]}], "4"], RowBox[{"IdentityMatrix", "[", RowBox[{"Length", "[", "system", "]"}], "]"}]}]}]}], "]"}]], "Input"], Cell["\<\ The level diagram showing resonant (co-rotating) optical couplings. \ \>", "Text"], Cell[BoxData[ RowBox[{"diagram", "=", RowBox[{"LevelDiagram", "[", RowBox[{"system", ",", RowBox[{ RowBox[{"Hrwa", "/.", RowBox[{"{", RowBox[{ RowBox[{ SubscriptBox["\[Delta]", "1"], "\[Rule]", RowBox[{"-", "\[Omega]1res"}]}], ",", RowBox[{ SubscriptBox["\[Delta]", "2"], "\[Rule]", RowBox[{"-", "\[Omega]2res"}]}], ",", RowBox[{ SubscriptBox["\[Delta]", "3"], "\[Rule]", RowBox[{"-", "\[Omega]3res"}]}]}], "}"}]}], "/.", RowBox[{"{", RowBox[{ RowBox[{ RowBox[{"Energy", "[", "0", "]"}], "\[Rule]", "0"}], ",", RowBox[{ RowBox[{"Energy", "[", "1", "]"}], "\[Rule]", "2"}], ",", RowBox[{ RowBox[{"Energy", "[", "2", "]"}], "\[Rule]", "3"}], ",", RowBox[{ RowBox[{"HyperfineA", "[", "_", "]"}], "\[Rule]", ".2"}], ",", RowBox[{ RowBox[{"HyperfineB", "[", "_", "]"}], "\[Rule]", ".1"}], ",", RowBox[{"\[CapitalOmega]L", "\[Rule]", ".1"}]}], "}"}]}], ",", RowBox[{"ParityOffset", "\[Rule]", "False"}], ",", RowBox[{"Epilog", "\[Rule]", RowBox[{"{", RowBox[{ RowBox[{"Text", "[", RowBox[{ RowBox[{"Style", "[", RowBox[{ "\"\<\!\(\*SuperscriptBox[\(\[VeryThinSpace]\), \ \(2\)]\)\!\(\*SubscriptBox[\(S\), \(1/2\)]\)\>\"", ",", RowBox[{"SingleLetterItalics", "\[Rule]", "False"}]}], "]"}], ",", RowBox[{"{", RowBox[{ RowBox[{"-", "2.9"}], ",", "0"}], "}"}]}], "]"}], ",", RowBox[{"Text", "[", RowBox[{ RowBox[{"Style", "[", RowBox[{ "\"\<\!\(\*SuperscriptBox[\(\[VeryThinSpace]\), \ \(2\)]\)\!\(\*SubscriptBox[\(P\), \(1/2\)]\)\>\"", ",", RowBox[{"SingleLetterItalics", "\[Rule]", "False"}]}], "]"}], ",", RowBox[{"{", RowBox[{ RowBox[{"-", "2.9"}], ",", "1.6"}], "}"}]}], "]"}], ",", RowBox[{"Text", "[", RowBox[{ RowBox[{"Style", "[", RowBox[{ "\"\<\!\(\*SuperscriptBox[\(\[VeryThinSpace]\), \ \(2\)]\)\!\(\*SubscriptBox[\(P\), \(3/2\)]\)\>\"", ",", RowBox[{"SingleLetterItalics", "\[Rule]", "False"}]}], "]"}], ",", RowBox[{"{", RowBox[{ RowBox[{"-", "2.9"}], ",", "2.7"}], "}"}]}], "]"}]}], "}"}]}], ",", RowBox[{"ImageSize", "\[Rule]", RowBox[{"3.5", " ", "72"}]}], ",", RowBox[{"ImagePadding", "\[Rule]", RowBox[{"{", RowBox[{ RowBox[{"{", RowBox[{"35", ",", "0"}], "}"}], ",", RowBox[{"{", RowBox[{"10", ",", "10"}], "}"}]}], "}"}]}]}], "]"}]}]], "Input", CellID->102096636], Cell[TextData[{ Cell[BoxData[ ButtonBox["IntrinsicRelaxation", BaseStyle->"Link", ButtonData->"paclet:AtomicDensityMatrix/ref/IntrinsicRelaxation"]]], " and ", Cell[BoxData[ ButtonBox["TransitRelaxation", BaseStyle->"Link", ButtonData->"paclet:AtomicDensityMatrix/ref/TransitRelaxation"]]], " supply the matrices describing relaxation due to spontaneous decay and \ atomic transit, respectively. ", Cell[BoxData[ FormBox[ SubscriptBox["\[Gamma]", "t"], TraditionalForm]]], " is the transit rate." }], "Text", CellID->610306692], Cell[BoxData[ RowBox[{"MatrixForm", "[", RowBox[{"relax", "=", RowBox[{ RowBox[{"IntrinsicRelaxation", "[", "system", "]"}], "+", RowBox[{"TransitRelaxation", "[", RowBox[{"system", ",", "\[Gamma]t"}], "]"}]}]}], "]"}]], "Input", CellID->645617687], Cell[TextData[{ Cell[BoxData[ ButtonBox["OpticalRepopulation", BaseStyle->"Link", ButtonData->"paclet:AtomicDensityMatrix/ref/OpticalRepopulation"]]], " and ", Cell[BoxData[ ButtonBox["TransitRepopulation", BaseStyle->"Link", ButtonData->"paclet:AtomicDensityMatrix/ref/TransitRepopulation"]]], " supply the matrices describing repopulation of the ground state due to \ spontaneous decay and atomic transit." }], "Text", CellID->854192725], Cell[BoxData[ RowBox[{"MatrixForm", "[", RowBox[{"repop", "=", RowBox[{ RowBox[{"OpticalRepopulation", "[", "system", "]"}], "+", RowBox[{"TransitRepopulation", "[", RowBox[{"system", ",", "\[Gamma]t"}], "]"}]}]}], "]"}]], "Input", CellID->465762594], Cell["Here are the evolution equations.", "Text", CellID->314466782], Cell[BoxData[ RowBox[{"TableForm", "[", RowBox[{"eqs0", "=", RowBox[{"LiouvilleEquation", "[", RowBox[{"system", ",", "Hrwa", ",", "relax", ",", "repop"}], "]"}]}], "]"}]], "Input"], Cell["\<\ Here we pull useful numerical atomic data for Rb out of the database. These \ numbers are in omega units, so that the unit \[OpenCurlyDoubleQuote]Hertz\ \[CloseCurlyDoubleQuote] actually corresponds to rad/s.\ \>", "Text"], Cell[BoxData[ RowBox[{"atomicdata", "=", RowBox[{"Join", "[", RowBox[{ RowBox[{"AtomicData", "[", RowBox[{ RowBox[{"{", RowBox[{"\"\\"", ",", "87", ",", RowBox[{"{", RowBox[{"\"\\"", ",", RowBox[{"{", RowBox[{"5", ",", "\"\\""}], "}"}]}], "}"}], ",", RowBox[{"{", RowBox[{"2", ",", "\"\\"", ",", FractionBox["1", "2"]}], "}"}]}], "}"}], ",", RowBox[{"{", "HyperfineA", "}"}], ",", RowBox[{"Label", "\[Rule]", "0"}]}], "]"}], ",", RowBox[{"AtomicData", "[", RowBox[{ RowBox[{"{", RowBox[{"\"\\"", ",", "87", ",", RowBox[{"{", RowBox[{"\"\\"", ",", RowBox[{"{", RowBox[{"5", ",", "\"\\""}], "}"}]}], "}"}], ",", RowBox[{"{", RowBox[{"2", ",", "\"\\"", ",", FractionBox["1", "2"]}], "}"}]}], "}"}], ",", RowBox[{"{", RowBox[{"HyperfineA", ",", "NaturalWidth"}], "}"}], ",", RowBox[{"Label", "\[Rule]", "1"}]}], "]"}], ",", RowBox[{"AtomicData", "[", RowBox[{ RowBox[{"{", RowBox[{"\"\\"", ",", "87", ",", RowBox[{"{", RowBox[{"\"\\"", ",", RowBox[{"{", RowBox[{"5", ",", "\"\\""}], "}"}]}], "}"}], ",", RowBox[{"{", RowBox[{"2", ",", "\"\\"", ",", FractionBox["3", "2"]}], "}"}]}], "}"}], ",", RowBox[{"{", RowBox[{"HyperfineA", ",", "HyperfineB", ",", "NaturalWidth"}], "}"}], ",", RowBox[{"Label", "\[Rule]", "2"}]}], "]"}], ",", RowBox[{"{", RowBox[{ RowBox[{ RowBox[{"Wavelength", "[", "1", "]"}], "->", RowBox[{"1", "/", RowBox[{"Wavenumber", "[", RowBox[{"{", RowBox[{"\"\\"", ",", "87", ",", RowBox[{"{", RowBox[{"\"\\"", ",", RowBox[{"{", RowBox[{"5", ",", "\"\\""}], "}"}]}], "}"}], ",", RowBox[{"{", RowBox[{"2", ",", "\"\\"", ",", FractionBox["1", "2"]}], "}"}]}], "}"}], "]"}]}]}], ",", RowBox[{ RowBox[{"Wavelength", "[", "2", "]"}], "->", RowBox[{"1", "/", RowBox[{"Wavenumber", "[", RowBox[{"{", RowBox[{"\"\\"", ",", "87", ",", RowBox[{"{", RowBox[{"\"\\"", ",", RowBox[{"{", RowBox[{"5", ",", "\"\\""}], "}"}]}], "}"}], ",", RowBox[{"{", RowBox[{"2", ",", "\"\\"", ",", FractionBox["3", "2"]}], "}"}]}], "}"}], "]"}]}]}]}], "}"}]}], "]"}]}]], "Input"], Cell[TextData[{ "Here we find DM elements that are always identically zero, so we can remove \ them from the evolution equations. We put in sample values for all of the \ parameters, use the values for the atomic data from above, and set time \ derivatives to zero with ", Cell[BoxData[ RowBox[{ RowBox[{ SuperscriptBox[ SubscriptBox["\[Rho]", RowBox[{"s1_", ",", "s2_"}]], "\[Prime]", MultilineFunction->None], "[", "t", "]"}], "\[Rule]", "0"}]]], " to find the steady state. " }], "Text"], Cell[BoxData[ RowBox[{ RowBox[{"sol", "=", RowBox[{ RowBox[{"NSolve", "[", RowBox[{ RowBox[{"Evaluate", "[", RowBox[{ RowBox[{ RowBox[{ RowBox[{ RowBox[{"eqs0", "/.", RowBox[{"{", RowBox[{ RowBox[{ SubscriptBox["\[Phi]", "_"], "\[Rule]", "1."}], ",", RowBox[{"\[Delta]k", "\[Rule]", "0"}], ",", RowBox[{"\[Delta]\[Omega]", "\[Rule]", "0"}], ",", RowBox[{ SubscriptBox["\[CapitalOmega]", "1"], "\[Rule]", SuperscriptBox["10", "11"]}], ",", RowBox[{ SubscriptBox["\[CapitalOmega]", "2"], "\[Rule]", RowBox[{"2.", " ", SuperscriptBox["10", "10"]}]}], ",", RowBox[{ SubscriptBox["\[CapitalOmega]", "_"], "\[Rule]", SuperscriptBox["10", "10"]}], ",", RowBox[{ SubscriptBox["\[Delta]", "_"], "\[Rule]", RowBox[{"1.", " ", SuperscriptBox["10", "7"]}]}], ",", RowBox[{"\[Gamma]t", "\[Rule]", "10."}], ",", RowBox[{"\[CapitalOmega]L", "\[Rule]", ".1"}]}], "}"}]}], "/.", "atomicdata"}], "/.", RowBox[{"Hertz", "\[Rule]", "1"}]}], "/.", RowBox[{"Mega", "\[Rule]", "1"}]}], "/.", RowBox[{ RowBox[{ SuperscriptBox[ SubscriptBox["\[Rho]", RowBox[{"s1_", ",", "s2_"}]], "\[Prime]", MultilineFunction->None], "[", "t", "]"}], "\[Rule]", "0"}]}], "]"}], ",", RowBox[{"DMVariables", "[", "system", "]"}]}], "]"}], "[", RowBox[{"[", "1", "]"}], "]"}]}], ";"}]], "Input"], Cell["\<\ Find the position in the solution list of all of the DM elements that are \ zero.\ \>", "Text"], Cell[BoxData[ RowBox[{"delpos", "=", RowBox[{"Position", "[", RowBox[{"sol", ",", RowBox[{"_", "\[Rule]", RowBox[{"0.", "+", RowBox[{"0.", " ", "\[ImaginaryI]"}]}]}]}], "]"}]}]], "Input"], Cell["\<\ Find the list of DM elements that correspond to these zero positions.\ \>", "Text"], Cell[BoxData[ RowBox[{ RowBox[{"delvars", "=", RowBox[{"Extract", "[", RowBox[{ RowBox[{"DMVariables", "[", "system", "]"}], ",", "delpos"}], "]"}]}], ";"}]], "Input"], Cell["Create a list of rules to set these DM elements to zero.", "Text"], Cell[BoxData[ RowBox[{ RowBox[{"delreps", "=", RowBox[{"(", RowBox[{ RowBox[{ RowBox[{"(", RowBox[{"#", "\[Rule]", "0"}], ")"}], "&"}], "/@", "delvars"}], ")"}]}], ";"}]], "Input"], Cell["\<\ Find the list of variables corresponding to nonzero DM elements.\ \>", "Text"], Cell[BoxData[ RowBox[{"vars", "=", RowBox[{"Delete", "[", RowBox[{ RowBox[{"DMVariables", "[", "system", "]"}], ",", "delpos"}], "]"}]}]], "Input"], Cell["\<\ Remove all of the zero elements from the evolution equations.\ \>", "Text"], Cell[BoxData[ RowBox[{"TableForm", "[", RowBox[{"eqs", "=", RowBox[{"Delete", "[", RowBox[{ RowBox[{"eqs0", "/.", "delreps"}], ",", "delpos"}], "]"}]}], "]"}]], "Input"], Cell["Initial conditions for the time-dependent case.", "Text"], Cell[BoxData[ RowBox[{"inits", "=", RowBox[{"Delete", "[", RowBox[{ RowBox[{"InitialConditions", "[", RowBox[{"system", ",", RowBox[{"TransitRepopulation", "[", RowBox[{"system", ",", "1"}], "]"}], ",", "0"}], "]"}], ",", "delpos"}], "]"}]}]], "Input"], Cell[TextData[{ "Find equations for the steady state by setting time derivatives to zero. We \ also substitute in numerical values for natural widths in units of ", Cell[BoxData[ FormBox[ SuperscriptBox["10", "6"], TraditionalForm]]], " rad/s." }], "Text"], Cell[BoxData[ RowBox[{"steadyeqs", "=", RowBox[{ RowBox[{ RowBox[{ RowBox[{ RowBox[{"eqs", "/.", RowBox[{ SubscriptBox["\[Phi]", "_"], "\[Rule]", "0"}]}], "/.", "atomicdata"}], "/.", RowBox[{"Hertz", "\[Rule]", "1"}]}], "/.", RowBox[{"Mega", "\[Rule]", "1"}]}], "/.", RowBox[{ RowBox[{ SuperscriptBox[ SubscriptBox["\[Rho]", RowBox[{"s1_", ",", "s2_"}]], "\[Prime]", MultilineFunction->None], "[", "t", "]"}], "\[Rule]", "0"}]}]}]], "Input"], Cell[TextData[{ "Here we find expressions for the fractional absorption (change in \ electric-field amplitude) and phase shift of the probe light per unit path \ length and unit atomic density, i.e., ", Cell[BoxData[ FormBox[ RowBox[{"{", RowBox[{ FractionBox["1", SubscriptBox["\[ScriptCapitalE]", "0"]], FractionBox[ SubscriptBox["d\[ScriptCapitalE]", "0"], RowBox[{ SubscriptBox["n", "0"], " ", "d\[ScriptL]"}]]}]}], TraditionalForm]]], ",", Cell[BoxData[ FormBox[ FractionBox["d\[Phi]", RowBox[{ SubscriptBox["n", "0"], " ", "d\[ScriptL]"}]], TraditionalForm]]], "}. These correspond to the imaginary and real parts of the index of \ refraction, respectively. Since the probe light does not interact with the D2 \ transition, we set DM elements involving the ", Cell[BoxData[ FormBox[ SubscriptBox["P", RowBox[{"3", "/", "2"}]], TraditionalForm]]], " state to zero. We also set the vanishing DM elements found above to zero. \ We find the expressions for arbitrary ellipticity and then set ellipticity to \ \[Pi]/4 after simplification in order to avoid an artificial divide by zero \ problem." }], "Text"], Cell[BoxData[{ RowBox[{ RowBox[{ RowBox[{"Observables", "[", RowBox[{"system", ",", RowBox[{"Energy", "[", "1", "]"}], ",", RowBox[{ SubscriptBox["\[CapitalOmega]", "2"], "/", RowBox[{"ReducedME", "[", RowBox[{"0", ",", RowBox[{"{", RowBox[{"Dipole", ",", "1"}], "}"}], ",", "1"}], "]"}]}], ",", RowBox[{"{", RowBox[{"0", ",", "\[Epsilon]"}], "}"}]}], "]"}], "[", RowBox[{"[", "1", "]"}], "]"}], ";"}], "\n", RowBox[{"obs0", "=", RowBox[{ RowBox[{"Simplify", "[", RowBox[{ RowBox[{"%", "/.", RowBox[{ RowBox[{ RowBox[{"DMElementPattern", "[", "]"}], "/;", RowBox[{ RowBox[{"Label2", "\[Equal]", "2"}], "||", RowBox[{"F1", "\[Equal]", "2"}]}]}], "\[Rule]", "0"}]}], "/.", "delreps"}], "]"}], "/.", RowBox[{"\[Epsilon]", "\[Rule]", RowBox[{"\[Pi]", "/", "4"}]}]}]}]}], "Input"], Cell["\<\ Calculate linear probe absorption with weak probe light on resonance and no \ pump to find reference level\ \>", "Text"], Cell[BoxData[{ RowBox[{ RowBox[{"params", "=", RowBox[{"{", RowBox[{ RowBox[{ SubscriptBox["\[CapitalOmega]", "1"], "\[Rule]", "0"}], ",", RowBox[{ SubscriptBox["\[CapitalOmega]", "2"], "\[Rule]", ".00001"}], ",", RowBox[{ SubscriptBox["\[CapitalOmega]", "3"], "\[Rule]", "0"}], ",", RowBox[{ SubscriptBox["\[CapitalOmega]", "4"], "\[Rule]", "0"}], ",", RowBox[{ SubscriptBox["\[Delta]", "2"], "\[Rule]", "0"}], ",", RowBox[{ SubscriptBox["\[Delta]", "1"], "\[Rule]", "0"}], ",", RowBox[{ SubscriptBox["\[Delta]", "3"], "\[Rule]", "0"}], ",", RowBox[{"\[Gamma]t", "\[Rule]", ".0001"}], ",", RowBox[{"\[CapitalOmega]L", "\[Rule]", "0"}]}], "}"}]}], ";"}], "\[IndentingNewLine]", RowBox[{"linabs", "=", RowBox[{ RowBox[{ RowBox[{"-", RowBox[{"obs0", "[", RowBox[{"[", "1", "]"}], "]"}]}], "/.", RowBox[{"Chop", "@", RowBox[{ RowBox[{"NSolve", "[", RowBox[{ RowBox[{"Evaluate", "[", RowBox[{"steadyeqs", "/.", "params"}], "]"}], ",", "vars"}], "]"}], "\[LeftDoubleBracket]", "1", "\[RightDoubleBracket]"}]}]}], "/.", "params"}]}]}], "Input"], Cell["\<\ Divide fractional absorption and phase shift by the linear absorption per \ unit length to find observables per absorption length.\ \>", "Text"], Cell[BoxData[ RowBox[{"obs", "=", RowBox[{ RowBox[{"{", RowBox[{ RowBox[{"-", "1"}], ",", "1"}], "}"}], RowBox[{"obs0", "/", "linabs"}]}]}]], "Input"] }, Open ]], Cell[CellGroupData[{ Cell["Plot steady-state results", "Section"], Cell[BoxData[ RowBox[{ RowBox[{"SetOptions", "[", RowBox[{"ListLinePlot", ",", RowBox[{"Frame", "\[Rule]", "True"}], ",", RowBox[{"Axes", "\[Rule]", "False"}], ",", RowBox[{"PlotRange", "\[Rule]", "All"}]}], "]"}], ";"}]], "Input"], Cell["\<\ Here we set parameters and plot the fractional absorption and phase shift of \ the probe light per absorption length, corresponding to imaginary and real \ parts of the index of refraction, respectively, as a function of probe \ detuning.\ \>", "Text"], Cell[BoxData[{ RowBox[{ RowBox[{"params", "=", RowBox[{"{", RowBox[{ RowBox[{ SubscriptBox["\[CapitalOmega]", "1"], "\[Rule]", "100."}], ",", RowBox[{ SubscriptBox["\[CapitalOmega]", "2"], "\[Rule]", ".01"}], ",", RowBox[{ SubscriptBox["\[CapitalOmega]", "3"], "\[Rule]", "10."}], ",", RowBox[{ SubscriptBox["\[CapitalOmega]", "4"], "\[Rule]", "0."}], ",", RowBox[{ SubscriptBox["\[Delta]", "2"], "\[Rule]", "d"}], ",", RowBox[{ SubscriptBox["\[Delta]", "1"], "\[Rule]", "0"}], ",", RowBox[{ SubscriptBox["\[Delta]", "3"], "\[Rule]", "0"}], ",", RowBox[{"\[Gamma]t", "\[Rule]", ".1"}], ",", RowBox[{"\[CapitalOmega]L", "\[Rule]", "0."}]}], "}"}]}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"steadyeqs1", "=", RowBox[{"steadyeqs", "/.", "params"}]}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"table", "=", RowBox[{"Table", "[", RowBox[{ RowBox[{ RowBox[{ RowBox[{"{", RowBox[{"d", ",", "obs"}], "}"}], "/.", RowBox[{"Chop", "@", RowBox[{ RowBox[{"NSolve", "[", RowBox[{"steadyeqs1", ",", "vars"}], "]"}], "\[LeftDoubleBracket]", "1", "\[RightDoubleBracket]"}]}]}], "/.", "params"}], ",", RowBox[{"{", RowBox[{"d", ",", RowBox[{"-", "200"}], ",", "200", ",", "2"}], "}"}]}], "]"}]}], ";"}], "\[IndentingNewLine]", RowBox[{"ListLinePlot", "[", RowBox[{"Transpose", "[", RowBox[{"{", RowBox[{ RowBox[{"table", "[", RowBox[{"[", RowBox[{"All", ",", "1"}], "]"}], "]"}], ",", RowBox[{"table", "[", RowBox[{"[", RowBox[{"All", ",", "2", ",", "1"}], "]"}], "]"}]}], "}"}], "]"}], "]"}], "\[IndentingNewLine]", RowBox[{"ListLinePlot", "[", RowBox[{"Transpose", "[", RowBox[{"{", RowBox[{ RowBox[{"table", "[", RowBox[{"[", RowBox[{"All", ",", "1"}], "]"}], "]"}], ",", RowBox[{"table", "[", RowBox[{"[", RowBox[{"All", ",", "2", ",", "2"}], "]"}], "]"}]}], "}"}], "]"}], "]"}]}], "Input"], Cell["Plot atomic populations for a given set of parameters.", "Text"], Cell[BoxData[{ RowBox[{ RowBox[{"params", "=", RowBox[{"{", RowBox[{ RowBox[{ SubscriptBox["\[CapitalOmega]", "1"], "\[Rule]", "10"}], ",", RowBox[{ SubscriptBox["\[CapitalOmega]", "2"], "\[Rule]", ".01"}], ",", RowBox[{ SubscriptBox["\[CapitalOmega]", "3"], "\[Rule]", "1."}], ",", RowBox[{ SubscriptBox["\[CapitalOmega]", "4"], "\[Rule]", "0."}], ",", RowBox[{ SubscriptBox["\[Delta]", "2"], "\[Rule]", "0"}], ",", RowBox[{ SubscriptBox["\[Delta]", "1"], "\[Rule]", "0"}], ",", RowBox[{ SubscriptBox["\[Delta]", "3"], "\[Rule]", "0"}], ",", RowBox[{"\[Gamma]t", "\[Rule]", ".1"}], ",", RowBox[{"\[CapitalOmega]L", "\[Rule]", "0."}]}], "}"}]}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"steadyeqs1", "=", RowBox[{"steadyeqs", "/.", "params"}]}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"sol", "=", RowBox[{ RowBox[{"Chop", "@", RowBox[{ RowBox[{"NSolve", "[", RowBox[{"steadyeqs1", ",", "vars"}], "]"}], "\[LeftDoubleBracket]", "1", "\[RightDoubleBracket]"}]}], "/.", "params"}]}], ";"}], "\n", RowBox[{"LevelDiagram", "[", RowBox[{"system", ",", RowBox[{ RowBox[{ RowBox[{"Hrwa", "/.", RowBox[{"{", RowBox[{ RowBox[{ SubscriptBox["\[Delta]", "1"], "\[Rule]", RowBox[{"-", "\[Omega]1res"}]}], ",", RowBox[{ SubscriptBox["\[Delta]", "2"], "\[Rule]", RowBox[{"-", "\[Omega]2res"}]}], ",", RowBox[{ SubscriptBox["\[Delta]", "3"], "\[Rule]", RowBox[{"-", "\[Omega]3res"}]}]}], "}"}]}], "/.", "params"}], "/.", RowBox[{"{", RowBox[{ RowBox[{ RowBox[{"Energy", "[", "0", "]"}], "\[Rule]", "0"}], ",", RowBox[{ RowBox[{"Energy", "[", "1", "]"}], "\[Rule]", "2"}], ",", RowBox[{ RowBox[{"Energy", "[", "2", "]"}], "\[Rule]", "3"}], ",", RowBox[{ RowBox[{"HyperfineA", "[", "_", "]"}], "\[Rule]", ".2"}], ",", RowBox[{ RowBox[{"HyperfineB", "[", "_", "]"}], "\[Rule]", ".1"}], ",", RowBox[{"\[CapitalOmega]L", "\[Rule]", "0"}]}], "}"}]}], ",", RowBox[{ RowBox[{ RowBox[{"DensityMatrix", "[", "system", "]"}], "/.", "sol"}], "/.", RowBox[{ RowBox[{"DMElementPattern", "[", "]"}], "\[Rule]", "0"}]}], ",", RowBox[{"ParityOffset", "\[Rule]", "False"}], ",", RowBox[{"Epilog", "\[Rule]", RowBox[{"{", RowBox[{ RowBox[{"Text", "[", RowBox[{ RowBox[{"Style", "[", RowBox[{ "\"\<\!\(\*SuperscriptBox[\(\[VeryThinSpace]\), \ \(2\)]\)\!\(\*SubscriptBox[\(S\), \(1/2\)]\)\>\"", ",", RowBox[{"SingleLetterItalics", "\[Rule]", "False"}]}], "]"}], ",", RowBox[{"{", RowBox[{ RowBox[{"-", "2.9"}], ",", "0"}], "}"}]}], "]"}], ",", RowBox[{"Text", "[", RowBox[{ RowBox[{"Style", "[", RowBox[{ "\"\<\!\(\*SuperscriptBox[\(\[VeryThinSpace]\), \ \(2\)]\)\!\(\*SubscriptBox[\(P\), \(1/2\)]\)\>\"", ",", RowBox[{"SingleLetterItalics", "\[Rule]", "False"}]}], "]"}], ",", RowBox[{"{", RowBox[{ RowBox[{"-", "2.9"}], ",", "1.6"}], "}"}]}], "]"}], ",", RowBox[{"Text", "[", RowBox[{ RowBox[{"Style", "[", RowBox[{ "\"\<\!\(\*SuperscriptBox[\(\[VeryThinSpace]\), \ \(2\)]\)\!\(\*SubscriptBox[\(P\), \(3/2\)]\)\>\"", ",", RowBox[{"SingleLetterItalics", "\[Rule]", "False"}]}], "]"}], ",", RowBox[{"{", RowBox[{ RowBox[{"-", "2.9"}], ",", "2.7"}], "}"}]}], "]"}]}], "}"}]}], ",", RowBox[{"ImageSize", "\[Rule]", RowBox[{"3.5", " ", "72"}]}], ",", RowBox[{"ImagePadding", "\[Rule]", RowBox[{"{", RowBox[{ RowBox[{"{", RowBox[{"35", ",", "0"}], "}"}], ",", RowBox[{"{", RowBox[{"10", ",", "10"}], "}"}]}], "}"}]}], ",", RowBox[{"PopulationStyle", "\[Rule]", RowBox[{"PointSize", "[", ".018", "]"}]}]}], "]"}]}], "Input"] }, Open ]], Cell[CellGroupData[{ Cell["Eliminate conjugate DM variables", "Section"], Cell[BoxData[{ RowBox[{"delvars1", "=", RowBox[{"Intersection", "[", RowBox[{"vars", ",", RowBox[{"Cases", "[", RowBox[{"vars", ",", RowBox[{ RowBox[{ RowBox[{"DMElementPattern", "[", "]"}], "/;", RowBox[{"!", RowBox[{"OrderedQ", "[", RowBox[{"{", RowBox[{"State2", ",", "State1"}], "}"}], "]"}]}]}], "\[Rule]", RowBox[{ SubscriptBox["\[Rho]", RowBox[{"State2", ",", "State1"}]], "[", "t", "]"}]}]}], "]"}]}], "]"}]}], "\[IndentingNewLine]", RowBox[{"delpos1", "=", RowBox[{"Flatten", "[", RowBox[{ RowBox[{ RowBox[{ RowBox[{"Position", "[", RowBox[{"vars", ",", "#"}], "]"}], "&"}], "/@", "%"}], ",", "1"}], "]"}]}], "\[IndentingNewLine]", RowBox[{"vars1", "=", RowBox[{"Delete", "[", RowBox[{"vars", ",", "delpos1"}], "]"}]}]}], "Input"], Cell[BoxData[ RowBox[{"delreps1", "=", RowBox[{"delvars1", "/.", RowBox[{ RowBox[{ SubscriptBox["\[Rho]", RowBox[{"s1_", ",", "s2_"}]], "[", "t", "]"}], "\[Rule]", RowBox[{"(", RowBox[{ RowBox[{ SubscriptBox["\[Rho]", RowBox[{"s1", ",", "s2"}]], "[", "t", "]"}], "\[Rule]", RowBox[{"Conjugate", "[", RowBox[{ SubscriptBox["\[Rho]", RowBox[{"s2", ",", "s1"}]], "[", "t", "]"}], "]"}]}], ")"}]}]}]}]], "Input"], Cell[BoxData[ RowBox[{"eqs1", "=", RowBox[{ RowBox[{"Delete", "[", RowBox[{"eqs", ",", "delpos1"}], "]"}], "/.", "delreps1"}]}]], "Input"], Cell[BoxData[ RowBox[{"inits1", "=", RowBox[{"Delete", "[", RowBox[{"inits", ",", "delpos1"}], "]"}]}]], "Input"] }, Open ]], Cell[CellGroupData[{ Cell["Convert equations to XMDS form", "Section"], Cell[BoxData[{ RowBox[{ RowBox[{ RowBox[{ RowBox[{ RowBox[{ RowBox[{"Join", "[", RowBox[{"atomicdata", ",", RowBox[{"{", RowBox[{ RowBox[{"rt6", "==", RowBox[{"N", "[", SqrtBox["6"], "]"}]}], ",", RowBox[{"rt3", "==", RowBox[{"N", "[", SqrtBox["3"], "]"}]}], ",", RowBox[{"rt2", "==", RowBox[{"N", "[", SqrtBox["2"], "]"}]}]}], "}"}]}], "]"}], "/.", RowBox[{"Mega", "\[Rule]", SuperscriptBox["10", "6"]}]}], "/.", RowBox[{"Hertz", "\[Rule]", "1"}]}], "/.", RowBox[{"Centimeter", "\[Rule]", "1"}]}], "/.", RowBox[{"Rule", "\[Rule]", "Equal"}]}], "/.", RowBox[{"{", RowBox[{ RowBox[{"NaturalWidth", "\[Rule]", "g"}], ",", RowBox[{"HyperfineA", "\[Rule]", "ha"}], ",", RowBox[{"HyperfineB", "\[Rule]", "hb"}], ",", RowBox[{"Wavelength", "\[Rule]", "lambda"}]}], "}"}]}], "\[IndentingNewLine]", RowBox[{"StringReplace", "[", RowBox[{ RowBox[{"StringJoin", "[", RowBox[{ RowBox[{ RowBox[{ RowBox[{"ToString", "@", RowBox[{"CForm", "[", "#", "]"}]}], "<>", "\"\<;\\n\>\""}], "&"}], "/@", "%"}], "]"}], ",", RowBox[{"{", RowBox[{ RowBox[{"\"\<==\>\"", "\[Rule]", "\"\<=\>\""}], ",", RowBox[{ RowBox[{"Shortest", "[", RowBox[{ "\"\\"", "~~", "a__", "~~", "\"\<,\>\"", "~~", "b__", "~~", "\"\<)\>\""}], "]"}], ":>", RowBox[{"\"\\"", "<>", "a", "<>", "b", "<>", "\"\<_dt\>\""}]}], ",", RowBox[{ RowBox[{"Shortest", "[", RowBox[{ "\"\\"", "~~", "a__", "~~", "\"\<,\>\"", "~~", "b__", "~~", "\"\<)\>\""}], "]"}], ":>", RowBox[{"\"\\"", "<>", "a", "<>", "b"}]}], ",", RowBox[{ RowBox[{"\"\\"", "~~", "a_", "~~", "\"\<)\>\""}], "\[RuleDelayed]", RowBox[{"\"\\"", "<>", "a"}]}], ",", RowBox[{ RowBox[{"\"\\"", "~~", "a_", "~~", "\"\<)\>\""}], "\[RuleDelayed]", RowBox[{"\"\\"", "<>", "a"}]}], ",", RowBox[{ RowBox[{"\"\\"", "~~", "a_", "~~", "\"\<)\>\""}], "\[RuleDelayed]", RowBox[{"\"\\"", "<>", "a"}]}], ",", RowBox[{ RowBox[{"\"\\"", "~~", "a_", "~~", "\"\<)\>\""}], "\[RuleDelayed]", RowBox[{"\"\\"", "<>", "a"}]}], ",", RowBox[{ RowBox[{"\"\\"", "~~", "a_", "~~", "\"\<)\>\""}], "\[RuleDelayed]", RowBox[{"\"\\"", "<>", "a"}]}], ",", RowBox[{ RowBox[{"\"\\"", "~~", "a_", "~~", "\"\<)\>\""}], "\[RuleDelayed]", RowBox[{"\"\\"", "<>", "a"}]}], ",", RowBox[{ RowBox[{"\"\\"", "~~", "a_", "~~", "\"\<)\>\""}], "\[RuleDelayed]", RowBox[{"\"\\"", "<>", "a"}]}], ",", " ", RowBox[{ "\"\\"", "\[RuleDelayed]", " ", "\"\\""}]}], "}"}]}], "]"}]}], "Input", CellChangeTimes->{{3.557577613086402*^9, 3.557577705341272*^9}, { 3.557577737830722*^9, 3.557577743453672*^9}, {3.557577791685976*^9, 3.557577804713066*^9}, {3.55757785531134*^9, 3.557577869097969*^9}, { 3.557577909846743*^9, 3.557577914049493*^9}, {3.557577958346201*^9, 3.557577990666655*^9}}], Cell[BoxData[ RowBox[{"Export", "[", RowBox[{ RowBox[{"ToFileName", "[", RowBox[{ RowBox[{"NotebookDirectory", "[", "]"}], ",", "\"\\""}], "]"}], ",", "%"}], "]"}]], "Input"], Cell["Convert equations to C form", "Text"], Cell[BoxData[{ RowBox[{ RowBox[{ RowBox[{"{", RowBox[{ RowBox[{"Label", "[", "#", "]"}], ",", RowBox[{"F", "[", "#", "]"}], ",", RowBox[{"M", "[", "#", "]"}]}], "}"}], "&"}], "/@", "system"}], "\[IndentingNewLine]", RowBox[{"labelrep", "=", RowBox[{"Table", "[", RowBox[{ RowBox[{ RowBox[{"%", "[", RowBox[{"[", "i", "]"}], "]"}], "\[Rule]", RowBox[{"ToString", "@", RowBox[{"PaddedForm", "[", RowBox[{"i", ",", "2", ",", RowBox[{"NumberPadding", "\[Rule]", RowBox[{"{", RowBox[{"\"\<0\>\"", ",", "\"\<\>\""}], "}"}]}], ",", RowBox[{"NumberSigns", "\[Rule]", RowBox[{"{", RowBox[{"\"\<\>\"", ",", "\"\<\>\""}], "}"}]}]}], "]"}]}]}], ",", RowBox[{"{", RowBox[{"i", ",", RowBox[{"Length", "@", "%"}]}], "}"}]}], "]"}]}]}], "Input"], Cell[BoxData[{ RowBox[{ RowBox[{ RowBox[{ RowBox[{ RowBox[{ RowBox[{ RowBox[{ RowBox[{"Join", "[", RowBox[{ RowBox[{"{", RowBox[{ RowBox[{"rt6", "==", RowBox[{"N", "[", SqrtBox["6"], "]"}]}], ",", RowBox[{"rt3", "==", RowBox[{"N", "[", SqrtBox["3"], "]"}]}], ",", RowBox[{"rt2", "==", RowBox[{"N", "[", SqrtBox["2"], "]"}]}]}], "}"}], ",", "eqs1"}], "]"}], "/.", RowBox[{"{", RowBox[{ RowBox[{ SubscriptBox["\[Delta]", "s_"], "\[Rule]", RowBox[{ SubscriptBox["\[Delta]", "s"], "+", RowBox[{ SubscriptBox["k", "s"], "v"}]}]}], ",", RowBox[{"\[Delta]\[Omega]", "\[Rule]", RowBox[{"\[Delta]\[Omega]", "+", RowBox[{"\[Delta]k", " ", "v"}]}]}]}], "}"}]}], "/.", RowBox[{"{", RowBox[{ RowBox[{"NaturalWidth", "\[Rule]", "g"}], ",", RowBox[{"HyperfineA", "\[Rule]", "ha"}], ",", RowBox[{"HyperfineB", "\[Rule]", "hb"}], ",", RowBox[{ RowBox[{ SuperscriptBox[ SubscriptBox["\[Rho]", RowBox[{"s1_", ",", "s2_"}]], "\[Prime]", MultilineFunction->None], "[", "t", "]"}], "\[RuleDelayed]", RowBox[{"dr", "[", RowBox[{"s1", ",", "s2"}], "]"}]}], ",", RowBox[{ RowBox[{ SubscriptBox["\[Rho]", RowBox[{"s1_", ",", "s2_"}]], "[", "t", "]"}], "\[RuleDelayed]", RowBox[{"r", "[", RowBox[{"s1", ",", "s2"}], "]"}]}], ",", RowBox[{"\[Gamma]t", "\[Rule]", "gt"}], ",", RowBox[{"\[CapitalOmega]L", "\[Rule]", "WL"}], ",", RowBox[{ SubscriptBox["\[Delta]", "i_"], "\[Rule]", RowBox[{"d", "[", "i", "]"}]}], ",", RowBox[{ SubscriptBox["k", "i_"], "\[Rule]", RowBox[{"k", "[", "i", "]"}]}], ",", RowBox[{ SuperscriptBox["\[ExponentialE]", RowBox[{"a_.", "+", RowBox[{"\[ImaginaryI]", " ", SubscriptBox["\[Phi]", "i_"]}]}]], "\[Rule]", RowBox[{ SuperscriptBox["\[ExponentialE]", "a"], RowBox[{ RowBox[{"W", "[", "i", "]"}], "/", SubscriptBox["\[CapitalOmega]", "i"]}]}]}], ",", RowBox[{ SuperscriptBox["\[ExponentialE]", RowBox[{"a_.", "-", RowBox[{"\[ImaginaryI]", " ", SubscriptBox["\[Phi]", "i_"]}]}]], "\[Rule]", RowBox[{ SuperscriptBox["\[ExponentialE]", "a"], RowBox[{ RowBox[{"Wc", "[", "i", "]"}], "/", SubscriptBox["\[CapitalOmega]", "i"]}]}]}]}], "}"}]}], "/.", "labelrep"}], "/.", RowBox[{ RowBox[{"Power", "[", RowBox[{"E", ",", "b_"}], "]"}], "\[Rule]", RowBox[{"exp", "[", "b", "]"}]}]}], "/.", RowBox[{"Conjugate", "\[Rule]", "conj"}]}], "/.", RowBox[{ RowBox[{"Complex", "[", RowBox[{"0", ",", "a_"}], "]"}], "\[Rule]", RowBox[{"a", " ", "i"}]}]}], "/.", RowBox[{"{", RowBox[{ RowBox[{ RowBox[{"1", "/", SqrtBox["6"]}], "\[Rule]", RowBox[{"1", "/", "rt6"}]}], ",", RowBox[{ RowBox[{"1", "/", SqrtBox["3"]}], "\[Rule]", RowBox[{"1", "/", "rt3"}]}], ",", RowBox[{ RowBox[{"1", "/", SqrtBox["2"]}], "\[Rule]", RowBox[{"1", "/", "rt2"}]}]}], "}"}]}], "\[IndentingNewLine]", RowBox[{"StringReplace", "[", RowBox[{ RowBox[{"StringJoin", "[", RowBox[{ RowBox[{ RowBox[{ RowBox[{"ToString", "@", RowBox[{"CForm", "[", "#", "]"}]}], "<>", "\"\<;\\n\>\""}], "&"}], "/@", "%"}], "]"}], ",", RowBox[{"{", RowBox[{ RowBox[{"\"\<==\>\"", "\[Rule]", "\"\<=\>\""}], ",", RowBox[{ RowBox[{"Shortest", "[", RowBox[{ "\"\\"", "~~", "a__", "~~", "\"\<\\\",\\\"\>\"", "~~", "b__", "~~", "\"\<\\\")\>\""}], "]"}], ":>", RowBox[{"\"\\"", "<>", "a", "<>", "b", "<>", "\"\<_dt\>\""}]}], ",", RowBox[{ RowBox[{"Shortest", "[", RowBox[{ "\"\\"", "~~", "a__", "~~", "\"\<\\\",\\\"\>\"", "~~", "b__", "~~", "\"\<\\\")\>\""}], "]"}], ":>", RowBox[{"\"\\"", "<>", "a", "<>", "b"}]}], ",", RowBox[{ RowBox[{"\"\\"", "~~", "a_", "~~", "\"\<)\>\""}], "\[RuleDelayed]", RowBox[{"\"\\"", "<>", "a"}]}], ",", RowBox[{ RowBox[{"\"\\"", "~~", "a_", "~~", "\"\<)\>\""}], "\[RuleDelayed]", RowBox[{"\"\\"", "<>", "a"}]}], ",", RowBox[{ RowBox[{"\"\\"", "~~", "a_", "~~", "\"\<)\>\""}], "\[RuleDelayed]", RowBox[{"\"\\"", "<>", "a"}]}], ",", RowBox[{ RowBox[{"\"\\"", "~~", "a_", "~~", "\"\<)\>\""}], "\[RuleDelayed]", RowBox[{"\"\\"", "<>", "a", "<>", "\"\\""}]}], ",", RowBox[{ RowBox[{"\"\\"", "~~", "a_", "~~", "\"\<)\>\""}], "\[RuleDelayed]", RowBox[{"\"\\"", "<>", "a", "<>", "\"\\""}]}], ",", RowBox[{ RowBox[{"\"\\"", "~~", "a_", "~~", "\"\<)\>\""}], "\[RuleDelayed]", RowBox[{"\"\\"", "<>", "a"}]}], ",", RowBox[{ RowBox[{"\"\\"", "~~", "a_", "~~", "\"\<)\>\""}], "\[RuleDelayed]", RowBox[{"\"\\"", "<>", "a"}]}]}], "}"}]}], "]"}]}], "Input"], Cell[BoxData[ RowBox[{"Export", "[", RowBox[{ RowBox[{"ToFileName", "[", RowBox[{ RowBox[{"NotebookDirectory", "[", "]"}], ",", "\"\\""}], "]"}], ",", "%"}], "]"}]], "Input"], Cell["Convert initial conditions to c form.", "Text"], Cell[BoxData[{ RowBox[{ RowBox[{ RowBox[{ RowBox[{ RowBox[{ RowBox[{"inits1", "/.", RowBox[{"{", RowBox[{ RowBox[{ SubscriptBox["\[Delta]", "s_"], "\[Rule]", RowBox[{ SubscriptBox["\[Delta]", "s"], "+", RowBox[{ SubscriptBox["k", "s"], "v"}]}]}], ",", RowBox[{"\[Delta]\[Omega]", "\[Rule]", RowBox[{"\[Delta]\[Omega]", "+", RowBox[{"\[Delta]k", " ", "v"}]}]}]}], "}"}]}], "/.", RowBox[{"{", RowBox[{ RowBox[{"NaturalWidth", "\[Rule]", "g"}], ",", RowBox[{"HyperfineA", "\[Rule]", "ha"}], ",", RowBox[{"HyperfineB", "\[Rule]", "hb"}], ",", RowBox[{ RowBox[{ SuperscriptBox[ SubscriptBox["\[Rho]", RowBox[{"s1_", ",", "s2_"}]], "\[Prime]", MultilineFunction->None], "[", "t_", "]"}], "\[RuleDelayed]", RowBox[{"dr", "[", RowBox[{"s1", ",", "s2"}], "]"}]}], ",", RowBox[{ RowBox[{ SubscriptBox["\[Rho]", RowBox[{"s1_", ",", "s2_"}]], "[", "t_", "]"}], "\[RuleDelayed]", RowBox[{"r", "[", RowBox[{"s1", ",", "s2"}], "]"}]}], ",", RowBox[{"\[Gamma]t", "\[Rule]", "gt"}], ",", RowBox[{ SubscriptBox["\[Delta]", "i_"], "\[Rule]", RowBox[{"d", "[", "i", "]"}]}], ",", RowBox[{ SubscriptBox["k", "i_"], "\[Rule]", RowBox[{"k", "[", "i", "]"}]}], ",", RowBox[{ SuperscriptBox["\[ExponentialE]", RowBox[{"a_.", "+", RowBox[{"\[ImaginaryI]", " ", SubscriptBox["\[Phi]", "i_"]}]}]], "\[Rule]", RowBox[{ SuperscriptBox["\[ExponentialE]", "a"], RowBox[{ RowBox[{"W", "[", "i", "]"}], "/", SubscriptBox["\[CapitalOmega]", "i"]}]}]}], ",", RowBox[{ SuperscriptBox["\[ExponentialE]", RowBox[{"a_.", "-", RowBox[{"\[ImaginaryI]", " ", SubscriptBox["\[Phi]", "i_"]}]}]], "\[Rule]", RowBox[{ SuperscriptBox["\[ExponentialE]", "a"], RowBox[{ RowBox[{"Wc", "[", "i", "]"}], "/", SubscriptBox["\[CapitalOmega]", "i"]}]}]}]}], "}"}]}], "/.", "labelrep"}], "/.", RowBox[{ RowBox[{"Power", "[", RowBox[{"E", ",", "b_"}], "]"}], "\[Rule]", RowBox[{"exp", "[", "b", "]"}]}]}], "/.", RowBox[{ RowBox[{"Complex", "[", RowBox[{"0", ",", "a_"}], "]"}], "\[Rule]", RowBox[{"a", " ", "i"}]}]}], "/.", RowBox[{"{", RowBox[{ RowBox[{ RowBox[{"1", "/", SqrtBox["6"]}], "\[Rule]", RowBox[{"1", "/", "rt6"}]}], ",", RowBox[{ RowBox[{"1", "/", SqrtBox["3"]}], "\[Rule]", RowBox[{"1", "/", "rt3"}]}], ",", RowBox[{ RowBox[{"1", "/", SqrtBox["2"]}], "\[Rule]", RowBox[{"1", "/", "rt2"}]}]}], "}"}]}], "\[IndentingNewLine]", RowBox[{"StringReplace", "[", RowBox[{ RowBox[{"StringJoin", "[", RowBox[{ RowBox[{ RowBox[{ RowBox[{"ToString", "@", RowBox[{"CForm", "[", "#", "]"}]}], "<>", "\"\<;\\n\>\""}], "&"}], "/@", "%"}], "]"}], ",", RowBox[{"{", RowBox[{ RowBox[{"\"\<==\>\"", "\[Rule]", "\"\<=\>\""}], ",", RowBox[{ RowBox[{"Shortest", "[", RowBox[{ "\"\\"", "~~", "a__", "~~", "\"\<\\\",\\\"\>\"", "~~", "b__", "~~", "\"\<\\\")\>\""}], "]"}], ":>", RowBox[{"\"\\"", "<>", "a", "<>", "b", "<>", "\"\<_dt\>\""}]}], ",", RowBox[{ RowBox[{"Shortest", "[", RowBox[{ "\"\\"", "~~", "a__", "~~", "\"\<\\\",\\\"\>\"", "~~", "b__", "~~", "\"\<\\\")\>\""}], "]"}], ":>", RowBox[{"\"\\"", "<>", "a", "<>", "b"}]}], ",", RowBox[{ RowBox[{"\"\\"", "~~", "a_", "~~", "\"\<)\>\""}], "\[RuleDelayed]", RowBox[{"\"\\"", "<>", "a"}]}], ",", RowBox[{ RowBox[{"\"\\"", "~~", "a_", "~~", "\"\<)\>\""}], "\[RuleDelayed]", RowBox[{"\"\\"", "<>", "a"}]}], ",", RowBox[{ RowBox[{"\"\\"", "~~", "a_", "~~", "\"\<)\>\""}], "\[RuleDelayed]", RowBox[{"\"\\"", "<>", "a"}]}], ",", RowBox[{ RowBox[{"\"\\"", "~~", "a_", "~~", "\"\<)\>\""}], "\[RuleDelayed]", RowBox[{"\"\\"", "<>", "a"}]}], ",", RowBox[{ RowBox[{"\"\\"", "~~", "a_", "~~", "\"\<)\>\""}], "\[RuleDelayed]", RowBox[{"\"\\"", "<>", "a"}]}], ",", RowBox[{ RowBox[{"\"\\"", "~~", "a_", "~~", "\"\<)\>\""}], "\[RuleDelayed]", RowBox[{"\"\\"", "<>", "a"}]}], ",", RowBox[{ RowBox[{"\"\\"", "~~", "a_", "~~", "\"\<)\>\""}], "\[RuleDelayed]", RowBox[{"\"\\"", "<>", "a"}]}]}], "}"}]}], "]"}]}], "Input", CellChangeTimes->{{3.557577072451949*^9, 3.557577074196031*^9}}], Cell[BoxData[ RowBox[{"Export", "[", RowBox[{ RowBox[{"ToFileName", "[", RowBox[{ RowBox[{"NotebookDirectory", "[", "]"}], ",", "\"\\""}], "]"}], ",", "%"}], "]"}]], "Input"] }, Open ]], Cell[CellGroupData[{ Cell["Find propagation equations", "Section"], Cell["Here we find the propagation equations.", "Text"], Cell["\<\ First, for each field 1--4, extract the list of sublevels that the field is \ assumed to interact with.\ \>", "Text"], Cell[BoxData[ RowBox[{"transitioncriteria", "=", RowBox[{"{", RowBox[{ RowBox[{ RowBox[{ RowBox[{"Label", "\[Equal]", "0"}], "&&", RowBox[{"F", "\[Equal]", "2"}]}], "||", RowBox[{ RowBox[{"Label", "\[Equal]", "1"}], "&&", RowBox[{"F", "\[Equal]", "1"}]}]}], ",", RowBox[{ RowBox[{ RowBox[{"Label", "\[Equal]", "0"}], "&&", RowBox[{"F", "\[Equal]", "1"}]}], "||", RowBox[{ RowBox[{"Label", "\[Equal]", "1"}], "&&", RowBox[{"F", "\[Equal]", "1"}]}]}], ",", RowBox[{ RowBox[{ RowBox[{"Label", "\[Equal]", "0"}], "&&", RowBox[{"F", "\[Equal]", "1"}]}], "||", RowBox[{ RowBox[{"Label", "\[Equal]", "2"}], "&&", RowBox[{"F", "\[Equal]", "2"}]}]}], ",", RowBox[{ RowBox[{ RowBox[{"Label", "\[Equal]", "0"}], "&&", RowBox[{"F", "\[Equal]", "2"}]}], "||", RowBox[{ RowBox[{"Label", "\[Equal]", "2"}], "&&", RowBox[{"F", "\[Equal]", "2"}]}]}]}], "}"}]}]], "Input"], Cell[BoxData[ RowBox[{"subsys", "=", RowBox[{"Table", "[", RowBox[{ RowBox[{"SelectStates", "[", RowBox[{"system", ",", RowBox[{"transitioncriteria", "[", RowBox[{"[", "k", "]"}], "]"}]}], "]"}], ",", RowBox[{"{", RowBox[{"k", ",", RowBox[{"Length", "[", "transitioncriteria", "]"}]}], "}"}]}], "]"}]}]], "Input"], Cell["\<\ The reduced matrix elements corresponding to the transition for each field.\ \>", "Text"], Cell[BoxData[ RowBox[{ RowBox[{"rmes", "=", RowBox[{"{", RowBox[{ RowBox[{"ReducedME", "[", RowBox[{"0", ",", RowBox[{"{", RowBox[{"Dipole", ",", "1"}], "}"}], ",", "1"}], "]"}], ",", " ", RowBox[{"ReducedME", "[", RowBox[{"0", ",", RowBox[{"{", RowBox[{"Dipole", ",", "1"}], "}"}], ",", "1"}], "]"}], ",", " ", RowBox[{"ReducedME", "[", RowBox[{"0", ",", RowBox[{"{", RowBox[{"Dipole", ",", "1"}], "}"}], ",", "2"}], "]"}], ",", " ", RowBox[{"ReducedME", "[", RowBox[{"0", ",", RowBox[{"{", RowBox[{"Dipole", ",", "1"}], "}"}], ",", "2"}], "]"}]}], "}"}]}], ";"}]], "Input"], Cell[TextData[{ "The first-order wave equation is given in terms of the rotating-frame \ density matrix by\n\t", Cell[BoxData[ FormBox[ RowBox[{ RowBox[{ FractionBox[ RowBox[{"\[PartialD]", SubscriptBox[ StyleBox["E", FontWeight->"Bold", FontSlant->"Plain"], "0"]}], RowBox[{"\[PartialD]", "z"}]], "+", RowBox[{ FractionBox["1", "c"], FractionBox[ RowBox[{"\[PartialD]", SubscriptBox[ StyleBox["E", FontWeight->"Bold", FontSlant->"Plain"], "0"]}], RowBox[{"\[PartialD]", "t"}]]}]}], "=", RowBox[{ FractionBox[ RowBox[{"4", "\[Pi]", " ", "\[ImaginaryI]", " ", "\[Omega]", " ", "n"}], "c"], RowBox[{ UnderscriptBox["\[Sum]", RowBox[{"m", "<", "n"}]], RowBox[{ SubscriptBox["\[Rho]", RowBox[{"n", ",", "m"}]], SubscriptBox[ StyleBox["d", FontWeight->"Bold", FontSlant->"Plain"], RowBox[{"m", ",", "n"}]], " "}]}]}]}], TraditionalForm]]] }], "Text"], Cell[TextData[{ "First we calculate the right-hand side. Here are the covariant spherical \ components of ", Cell[BoxData[ FormBox[ RowBox[{ UnderscriptBox["\[Sum]", RowBox[{"m", "<", "n"}]], RowBox[{ SubscriptBox["\[Rho]", RowBox[{"n", ",", "m"}]], SubscriptBox[ StyleBox["d", FontWeight->"Bold", FontSlant->"Plain"], RowBox[{"m", ",", "n"}]]}]}], TraditionalForm]]], " for each of the four transitions, in terms of the reduced dipole matrix \ element ||d||." }], "Text"], Cell[BoxData[ RowBox[{"polarizationcomponents1", "=", RowBox[{ RowBox[{"Table", "[", RowBox[{ RowBox[{"Sum", "[", RowBox[{ RowBox[{ RowBox[{ RowBox[{"DensityMatrix", "[", RowBox[{ "subsys", "\[LeftDoubleBracket]", "k", "\[RightDoubleBracket]"}], "]"}], "\[LeftDoubleBracket]", RowBox[{"n", ",", "m"}], "\[RightDoubleBracket]"}], RowBox[{"WignerEckart", "[", RowBox[{ RowBox[{ RowBox[{ "subsys", "\[LeftDoubleBracket]", "k", "\[RightDoubleBracket]"}], "\[LeftDoubleBracket]", "m", "\[RightDoubleBracket]"}], ",", RowBox[{"{", RowBox[{"Dipole", ",", "1"}], "}"}], ",", RowBox[{ RowBox[{ "subsys", "\[LeftDoubleBracket]", "k", "\[RightDoubleBracket]"}], "\[LeftDoubleBracket]", "n", "\[RightDoubleBracket]"}]}], "]"}]}], ",", RowBox[{"{", RowBox[{"n", ",", RowBox[{"Length", "[", RowBox[{ "subsys", "\[LeftDoubleBracket]", "k", "\[RightDoubleBracket]"}], "]"}]}], "}"}], ",", RowBox[{"{", RowBox[{"m", ",", "n"}], "}"}]}], "]"}], ",", RowBox[{"{", RowBox[{"k", ",", RowBox[{"Length", "[", "subsys", "]"}]}], "}"}]}], "]"}], "/.", "delreps"}]}]], "Input"], Cell[TextData[{ "Multiply through by ", Cell[BoxData[ FormBox[ RowBox[{ RowBox[{"4", "\[Pi]", " ", "\[ImaginaryI]", " ", "\[Omega]", " ", "n"}], " ", "||", "d", "||"}], TraditionalForm]]], ", and rewrite ", Cell[BoxData[ FormBox[ RowBox[{"||", "d", SuperscriptBox["||", "2"]}], TraditionalForm]]], " in terms of the natural width \[CapitalGamma] of the transition." }], "Text"], Cell[BoxData[ RowBox[{"polarizationcomponents2", "=", RowBox[{"Table", "[", RowBox[{ RowBox[{"Expand", "[", RowBox[{"ExpandDipoleRME", "[", RowBox[{ RowBox[{ "subsys", "\[LeftDoubleBracket]", "k", "\[RightDoubleBracket]"}], ",", RowBox[{ RowBox[{ RowBox[{"4", "\[Pi]", " ", "\[ImaginaryI]", " ", RowBox[{"Energy", "[", RowBox[{"Last", "[", RowBox[{ "rmes", "\[LeftDoubleBracket]", "k", "\[RightDoubleBracket]"}], "]"}], "]"}], " ", "n0", " ", RowBox[{ "rmes", "\[LeftDoubleBracket]", "k", "\[RightDoubleBracket]"}], " ", RowBox[{ "polarizationcomponents1", "\[LeftDoubleBracket]", "k", "\[RightDoubleBracket]"}]}], "/.", "delreps"}], "/.", "delreps1"}]}], "]"}], "]"}], ",", RowBox[{"{", RowBox[{"k", ",", RowBox[{"Length", "[", "subsys", "]"}]}], "}"}]}], "]"}]}]], "Input"], Cell[TextData[{ "Make the substitution ", Cell[BoxData[ FormBox[ RowBox[{ SubscriptBox["\[Eta]", "s"], "=", FractionBox[ RowBox[{"3", SubsuperscriptBox["\[Lambda]", "s", "2"], SubscriptBox["\[CapitalGamma]", "s"], "n"}], RowBox[{"4", "\[Pi]"}]]}], TraditionalForm]]], "." }], "Text"], Cell[BoxData[ RowBox[{"polarizationcomponents3", "=", RowBox[{ RowBox[{ RowBox[{"polarizationcomponents2", "/.", RowBox[{ RowBox[{"Energy", "[", "s_", "]"}], "\[Rule]", RowBox[{"2", RowBox[{"\[Pi]", "/", RowBox[{"\[Lambda]", "[", "s", "]"}]}]}]}]}], "/.", RowBox[{ RowBox[{"NaturalWidth", "[", "s_", "]"}], "\[Rule]", RowBox[{"4", "\[VeryThinSpace]", "\[Pi]", " ", RowBox[{ RowBox[{"\[Eta]", "[", "s", "]"}], "/", RowBox[{"(", RowBox[{"n0", " ", SuperscriptBox[ RowBox[{"\[Lambda]", "[", "s", "]"}], "2"], "3"}], ")"}]}]}]}]}], "//", "Simplify"}]}]], "Input"], Cell["\<\ Cartesian components of the amplitude of each of the 4 fields.\ \>", "Text"], Cell[BoxData[ RowBox[{"fieldcomponents1", "=", RowBox[{"Table", "[", RowBox[{ RowBox[{"field", "/.", RowBox[{"{", RowBox[{ RowBox[{ SubscriptBox["k", "_"], "\[Rule]", "0"}], ",", RowBox[{"t", "\[Rule]", "0"}], ",", RowBox[{ SubscriptBox["\[Phi]", "_"], "\[Rule]", "0"}], ",", RowBox[{ RowBox[{ SubscriptBox["\[CapitalOmega]", "s_"], "/;", RowBox[{"i", "\[NotEqual]", "s"}]}], "\[Rule]", "0"}]}], "}"}]}], ",", RowBox[{"{", RowBox[{"i", ",", "4"}], "}"}]}], "]"}]}]], "Input"], Cell["\<\ Convert to covariant spherical components and multiply by ||d||.\ \>", "Text"], Cell[BoxData[ RowBox[{"fieldcomponents2", "=", RowBox[{ RowBox[{ RowBox[{"ToCovariant", "/@", "fieldcomponents1"}], " ", "rmes"}], "/.", RowBox[{ SubscriptBox["\[CapitalOmega]", "s_"], "->", RowBox[{ SubscriptBox["\[CapitalOmega]", "s"], "[", "z", "]"}]}]}]}]], "Input"], Cell["\<\ Find propagation equations using the spherical components of the polarization \ and fields. Here we leave out the time-derivative term, which we will add in \ below.\ \>", "Text"], Cell[BoxData[ RowBox[{"propeqs", "=", RowBox[{ RowBox[{ RowBox[{"Solve", "[", RowBox[{ RowBox[{ RowBox[{"D", "[", RowBox[{"fieldcomponents2", ",", "z"}], "]"}], "==", RowBox[{"Simplify", "@", "polarizationcomponents3"}]}], ",", RowBox[{"Table", "[", RowBox[{ RowBox[{ SuperscriptBox[ SubscriptBox["\[CapitalOmega]", "i"], "\[Prime]", MultilineFunction->None], "[", "z", "]"}], ",", RowBox[{"{", RowBox[{"i", ",", "4"}], "}"}]}], "]"}]}], "]"}], "[", RowBox[{"[", "1", "]"}], "]"}], "/.", RowBox[{"Rule", "\[Rule]", "Equal"}]}]}]], "Input"], Cell["Convert propagation equations to c form.", "Text"], Cell[BoxData[{ RowBox[{ RowBox[{ RowBox[{ RowBox[{ RowBox[{ RowBox[{"propeqs", "/.", RowBox[{"{", RowBox[{ RowBox[{"\[Eta]", "\[Rule]", "eta"}], ",", RowBox[{"NaturalWidth", "\[Rule]", "g"}], ",", RowBox[{ RowBox[{ SubscriptBox["\[Rho]", RowBox[{"s1_", ",", "s2_"}]], "[", "t", "]"}], "\[RuleDelayed]", RowBox[{"r", "[", RowBox[{"s1", ",", "s2"}], "]"}]}], ",", RowBox[{ RowBox[{ SuperscriptBox[ SubscriptBox["\[CapitalOmega]", "i_"], "\[Prime]", MultilineFunction->None], "[", "z", "]"}], "\[Rule]", RowBox[{"dW", "[", "i", "]"}]}]}], "}"}]}], "/.", "labelrep"}], "/.", RowBox[{ RowBox[{"Complex", "[", RowBox[{"0", ",", "a_"}], "]"}], "\[Rule]", RowBox[{"a", " ", "i"}]}]}], "/.", RowBox[{"Conjugate", "\[Rule]", "conj"}]}], "/.", RowBox[{ RowBox[{"a", ":", RowBox[{"Power", "[", RowBox[{ RowBox[{"_", "?", "NumericQ"}], ",", "_"}], "]"}]}], "\[RuleDelayed]", RowBox[{"N", "[", "a", "]"}]}]}], "//", "Simplify"}], "\[IndentingNewLine]", RowBox[{"StringJoin", "[", RowBox[{"Table", "[", RowBox[{ RowBox[{ RowBox[{"ToString", "@", RowBox[{"CForm", "[", RowBox[{"%", "[", RowBox[{"[", "j", "]"}], "]"}], "]"}]}], "<>", "\"\< - Lt[E\>\"", "<>", RowBox[{"ToString", "@", "j"}], "<>", "\"\<];\\n\>\""}], ",", RowBox[{"{", RowBox[{"j", ",", "4"}], "}"}]}], "]"}], "]"}], "\[IndentingNewLine]", RowBox[{"StringReplace", "[", RowBox[{"%", ",", RowBox[{"{", RowBox[{ RowBox[{"\"\<==\>\"", "\[Rule]", "\"\<=\>\""}], ",", RowBox[{ RowBox[{"\"\\"", "~~", "a_", "~~", "\"\<)\>\""}], "\[RuleDelayed]", RowBox[{"\"\\"", "<>", "a"}]}], ",", RowBox[{ RowBox[{"Shortest", "[", RowBox[{ "\"\\"", "~~", "a__", "~~", "\"\<\\\",\\\"\>\"", "~~", "b__", "~~", "\"\<\\\")\>\""}], "]"}], ":>", RowBox[{"\"\\"", "<>", "a", "<>", "b", "<>", "\"\<_dt\>\""}]}], ",", RowBox[{ RowBox[{"Shortest", "[", RowBox[{ "\"\\"", "~~", "a__", "~~", "\"\<\\\",\\\"\>\"", "~~", "b__", "~~", "\"\<\\\")\>\""}], "]"}], ":>", RowBox[{"\"\\"", "<>", "a", "<>", "b"}]}], ",", RowBox[{ RowBox[{"\"\\"", "~~", "a_", "~~", "\"\<)\>\""}], "\[RuleDelayed]", RowBox[{"\"\\"", "<>", "a"}]}], ",", RowBox[{ RowBox[{"\"\\"", "~~", "a_", "~~", "\"\<)\>\""}], "\[RuleDelayed]", RowBox[{"\"\\"", "<>", "a"}]}], ",", RowBox[{ RowBox[{"\"\\"", "~~", "a_", "~~", "\"\<)\>\""}], "\[RuleDelayed]", RowBox[{"\"\\"", "<>", "a"}]}], ",", RowBox[{ RowBox[{"\"\\"", "~~", "a_", "~~", "\"\<)\>\""}], "\[RuleDelayed]", RowBox[{"\"\\"", "<>", "a", "<>", "\"\<_dz\>\""}]}], ",", RowBox[{ RowBox[{"\"\\"", "~~", "a_", "~~", "\"\<)\>\""}], "\[RuleDelayed]", RowBox[{"\"\\"", "<>", "a"}]}], ",", RowBox[{ RowBox[{"\"\\"", "~~", "a_", "~~", "\"\<)\>\""}], "\[RuleDelayed]", RowBox[{"\"\\"", "<>", "a"}]}]}], "}"}]}], "]"}], "\[IndentingNewLine]", RowBox[{"Export", "[", RowBox[{ RowBox[{"ToFileName", "[", RowBox[{ RowBox[{"NotebookDirectory", "[", "]"}], ",", "\"\\""}], "]"}], ",", "%"}], "]"}]}], "Input"] }, Open ]] }, WindowSize->{956, 980}, WindowMargins->{{Automatic, 0}, {26, Automatic}}, PrivateNotebookOptions->{"FileOutlineCache"->False}, ShowSelection->True, FrontEndVersion->"8.0 for Linux x86 (32-bit) (February 23, 2011)", StyleDefinitions->"Default.nb" ]