Introduction to Cell and Real Atoms

The graph-based framework that rydiqule is based on is incredibly flexible, but often flexibilty is not as useful as being able to perform common actions quickly. For example, the Sensor module of rydiqule allows for any value of any parameter to be supported. While this is great for flexibility, if you know the first two states of an atomic system are all m_j fine structrue states for the D2 line of Rubidium-87, it would be somewhat painful to manually specify the dipole moment and transition frequencies for every possible state pairing for the laser coupling the the two manifolds.

The Cell module of rydiqule is designed to address exactly this issue. It builds upon Sensor with additional functionality to automatically add quantities like quantum numbers, dipole moments, state energies, and decoherence rates. Indeed, if you are familiar with object-oriented programming in python, Cell in fact inherits Sensor as a subclass, meaning it has access to all the same methods, and has the same internal underlying structure. Importantly, Cell is specifically designed to model alkali atoms supported by the Alkali Rydberg Calculator project.

This notebook can be downloaded here.

0. What is a Cell?

As we alluded to above, a Cell is just a Sensor at the end of the day. The familiar functions like add_coupling and get_hamiltonian are still there, but there is a lot more under the hood.

This notebook will go over some of the basic ways to use Cell to solve real systems quickly and easily. This tutorial assumes you have gone through and understand the basics laid out in the Introduction_To_Rydiqule notebook. If any of the terminolgy in this notebook is unfamiliar, revisit that notebook for a demonstration of basic principles. This notebook will not revisit those basic priciples, it will primarily highlight the differences between Cell and Sensor.

This notebook will also assume a basic familiarity with rydberg atom atomic physics. While we do have some physics documentation, rydiqule is aimed primarily at physicists who already understand many of these concepts, and the docs are more for our implementation of the underlying physics.

For starters, we will again import rydiqule as rq, as well as the usual imports. We also import A_QState directly for reasons that will become clear momentarily.

import rydiqule as rq
from rydiqule import A_QState

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

1. Creating a Cell

1.1 A_QState and constructor basics

If you remember in Sensor, the states could be integers, strings, or tuples. Cell is a little bit more restrictive about what you can use for states. To describe states in Cell, rydiqule has introduced a new type using python’s namedtuple called A_QState, short for atomic quantum state. Note that this is a subclass of tuple, so no rules have changed; using this is a special case of the bases that can be defined in Sensor. This named tuple has 3 mandatory arguments, n, l, and j representing the first 3 quantum number of a rydberg quantum state. Additionally, there are 3 more optional arguments for m_j, f, and m_f. As you might expect, you cannot specify f or m_f if m_j is specified, or vice versa. Also, to specify m_f, f must also be specified. All the usual quantum number rules for rydberg atoms must be obeyed.

We can access an A_QState directly as rq.A_QState, but it can be easier to from rydiqule import A_QState to shorten, since you may need to call it often. Let’s start with a simple ground and excited state for the D1 line of Rubidium-85 just to see how the Cell constructor works. As a minimum, we need to specify a atom with a string, and pass a list of A_QStates.

g = A_QState(5, 0, 0.5)
e = A_QState(5, 1, 0.5)

Rb_Cell = rq.Cell("Rb85", [g, e])
print(Rb_Cell.states)
print(type(Rb_Cell.states[0]))

[(n=5, l=0, j=0.5), (n=5, l=1, j=0.5)]
<class 'rydiqule.atom_utils.A_QState'>

The atom string can be any atoms supported by the ARC package, and following the <atomic symbol><isotope number> as above. Consult arc to see what isotopes are supported.

For readablity, A_QState trims off the "A_QState" text and unused quantum numbers from the string output, but rest assured, it is still an A_QState, which we see in the printout. You can see that the states of the system are exactly what are in the list we passed to the constructor. With a very simple system defined, we can inspect the Cell.couplings graph to see what sorts of things are on the graph.

print(Rb_Cell.couplings.nodes(data=True))
print(Rb_Cell.couplings.edges(data=True))

[((n=5, l=0, j=0.5), {'energy': 0.0, 'gamma_lifetime': 0.0}), ((n=5, l=1, j=0.5), {'energy': 2369435883.882498, 'gamma_lifetime': 36.11450417508357})]
[((n=5, l=1, j=0.5), (n=5, l=0, j=0.5), {'gamma_transition': 36.11450417508357, 'label': '((5, 1, 0.5),(5, 0, 0.5))'})]

We can see that a couple of things have been added to the graph. Each node has been populated automatically with energy levels (with the ground state defined at 0 Mrad), and state lifetimes. Like basically everything else in rydiqule, these quantities are expressed in Mrad/s. In addition to data on the nodes (which is just for reference and does not affect solving directly), there is a decoherent transition added from the second state to the first, associated with natural state lifetime of the \(5P_{1/2}\) state before it decays back to the ground state. In general, rydiqule will add natural transition rates from higher to lower states. There are some caveats and details to how this is done that will be discussed later in this section. Just for fun, let’s demonstrate that even without anything else explicitly defined, we already have a decoherence_matrix.

print(Rb_Cell.decoherence_matrix())

[[ 0.          0.        ]
 [36.11450418  0.        ]]

1.2 Efficiently defining A_QStates

List shorthand

While it is nice that we can define all of these states with A_QState, it is easy to imagine a situation in which there are a lot more states in a system of interest that you want to write down (consider all possible hyperfine sublevels in a high angular momentum upper rydberg state). You could be defining the states for hours. As you might imagine, rydiqule has a built-in way to define all of these states at once. When we define our A_QState, we can define any of the quantum numbers from j onwards (so m_j, f, m_f) as a list, expanding the single specification out into a list of states automatically. Below we show this shorthand to add both the D1 and D2 excited states with one A_QState specification. This functionality is identical to using tuple states in a Sensor, and indeed is the motivation for adding the functionality in the first place.

g = A_QState(5, 0, 0.5)
e = A_QState(5, 1, [0.5, 1.5])

Rb_Cell_D12 = rq.Cell("Rb85", [g, e])
print(Rb_Cell_D12.states)
print(Rb_Cell_D12.decoherence_matrix())

[(n=5, l=0, j=0.5), (n=5, l=1, j=0.5), (n=5, l=1, j=1.5)]
[[ 0.          0.          0.        ]
 [36.11450418  0.          0.        ]
 [38.11316014  0.          0.        ]]

We can see that again we have created a Cell with decoherence values from the natural state lifetime already added, and once again we get the decoherence_matrix of the system (with no additional broadening) automatically.

“All” shorthand

Even providing things as lists can get cumbersome once sublevels start to get involved. For this reason, rydiqule supports using the string "all" for quantum numbers rather than just a list, which will automatically get states containing all allowed values of the specified quantum numbers. Let us again consider just the D2 line, but account for the m_j splitting of these levels.

g = A_QState(5, 0, 0.5, m_j="all")
e = A_QState(5, 1, 1.5, m_j="all")

D2_Cell_mj = rq.Cell("Rb85", [g,e])
print(D2_Cell_mj.states)
print(D2_Cell_mj.decoherence_matrix())

[(n=5, l=0, j=0.5, m_j=-0.5), (n=5, l=0, j=0.5, m_j=0.5), (n=5, l=1, j=1.5, m_j=-1.5), (n=5, l=1, j=1.5, m_j=-0.5), (n=5, l=1, j=1.5, m_j=0.5), (n=5, l=1, j=1.5, m_j=1.5)]
[[ 0.          0.          0.          0.          0.          0.        ]
 [ 0.          0.          0.          0.          0.          0.        ]
 [38.11316014  0.          0.          0.          0.          0.        ]
 [25.40877343 12.70438671  0.          0.          0.          0.        ]
 [12.70438671 25.40877343  0.          0.          0.          0.        ]
 [ 0.         38.11316014  0.          0.          0.          0.        ]]

We can see that we have 2 states in the ground manifold (corresponding to \(m_j=\pm 0.5\)), and 4 in the excited manifold (corresponding to \(m_j=\pm 0.5, \pm 1.5\)). Once again, we have sucessfully accounted for the transitions from excited into the ground state, including zeros for transitions which are not dipole-allowed.

We can take this "all" notion one step further by defining an entire hyperfine manifold. This time we will account for hyperfine splitting in both the ground and excited state. Note when computing hyperfine splitting, since we supplied the atom specification to the constructor, rydiqule already knows the nuclear magnetic moment of our atom (Rubidium-85) to be \(\frac{5}{2}\), and will calculate allowed f values accordingly

g = A_QState(5, 0, 0.5, f="all", m_f="all")
e = A_QState(5, 1, 1.5, f="all", m_f="all")

D2_Cell_hyperfine = rq.Cell("Rb85", [g,e])
print(D2_Cell_hyperfine.states)
print()
print(f"{len(D2_Cell_hyperfine.states)} states")

[(n=5, l=0, j=0.5, f=2.0, m_f=-2.0), (n=5, l=0, j=0.5, f=2.0, m_f=-1.0), (n=5, l=0, j=0.5, f=2.0, m_f=0.0), (n=5, l=0, j=0.5, f=2.0, m_f=1.0), (n=5, l=0, j=0.5, f=2.0, m_f=2.0), (n=5, l=0, j=0.5, f=3.0, m_f=-3.0), (n=5, l=0, j=0.5, f=3.0, m_f=-2.0), (n=5, l=0, j=0.5, f=3.0, m_f=-1.0), (n=5, l=0, j=0.5, f=3.0, m_f=0.0), (n=5, l=0, j=0.5, f=3.0, m_f=1.0), (n=5, l=0, j=0.5, f=3.0, m_f=2.0), (n=5, l=0, j=0.5, f=3.0, m_f=3.0), (n=5, l=1, j=1.5, f=1.0, m_f=-1.0), (n=5, l=1, j=1.5, f=1.0, m_f=0.0), (n=5, l=1, j=1.5, f=1.0, m_f=1.0), (n=5, l=1, j=1.5, f=2.0, m_f=-2.0), (n=5, l=1, j=1.5, f=2.0, m_f=-1.0), (n=5, l=1, j=1.5, f=2.0, m_f=0.0), (n=5, l=1, j=1.5, f=2.0, m_f=1.0), (n=5, l=1, j=1.5, f=2.0, m_f=2.0), (n=5, l=1, j=1.5, f=3.0, m_f=-3.0), (n=5, l=1, j=1.5, f=3.0, m_f=-2.0), (n=5, l=1, j=1.5, f=3.0, m_f=-1.0), (n=5, l=1, j=1.5, f=3.0, m_f=0.0), (n=5, l=1, j=1.5, f=3.0, m_f=1.0), (n=5, l=1, j=1.5, f=3.0, m_f=2.0), (n=5, l=1, j=1.5, f=3.0, m_f=3.0), (n=5, l=1, j=1.5, f=4.0, m_f=-4.0), (n=5, l=1, j=1.5, f=4.0, m_f=-3.0), (n=5, l=1, j=1.5, f=4.0, m_f=-2.0), (n=5, l=1, j=1.5, f=4.0, m_f=-1.0), (n=5, l=1, j=1.5, f=4.0, m_f=0.0), (n=5, l=1, j=1.5, f=4.0, m_f=1.0), (n=5, l=1, j=1.5, f=4.0, m_f=2.0), (n=5, l=1, j=1.5, f=4.0, m_f=3.0), (n=5, l=1, j=1.5, f=4.0, m_f=4.0)]

36 states

Obviously, we get a lot of states when we account for hyperfine splitting. While we do show how easy something like this, this is also somewhat of a warning. On the one hand, it is easy to account for tons of states with just a couple lines. On the other hand, a larger basis will take longer to solve (only in polynomial time, but still longer). This is not to discourage you from using the "all" feature, but to point out that you should only do so when you are looking to model physics only available when accounting for hyperfine. A couple of notes about some rules of hyperfine splitting:

1. We do not support mixing n,l,j states with hyperfine or fine states in the same Cell. If you want to use hyperfine, all states must be either fine split or hyperfine split.

2. A_QStates cannot include f without m_f. If you want hyperfine splitting and specify f, m_f must be a single value, list of allowed values, or "all". In any case, rydiqule will enforce that these quantum numbers are physically allowed.

Helper function shorthand

This is less crucial, but one more way rydiqule provides to specify manifolds of states in Cell is with a handful of utility functions that return lists of states. These can be useful if you either can’t be bothered to type the full A_QState out yourself or if you want to test a system with multiple different atoms and only change a single value in code. The relevant functions are as follows. In each case, n can be either an integer or a string of a particular atom as specified in the constructor (eg "Rb85"). splitting is one of [None, "fs", or "hfs"], with the default being None:

1. rq.ground_state(n, splitting=...), gets the \(l=0\), \(j=0.5\) state(s). This is what the in the “energy” value on the nodes seen previouslt are in reference to.

2. rq.D1_excited(n, splitting=...), gets the \(l=1\), \(j=0.5\) state(s).

3. rq.D2_excited(n, splitting=...), gets the \(l=1\), \(j=1.5\) state(s).

4. rq.D1_states(n, splitting=..., g_splitting=..., e_splitting=..) just calls 1 and 2.splitting overrrides ground and excited splitting if present.

5. rq.D2_states(n, splitting=..., g_splitting=..., e_splitting=..) just calls 1 and 2. splitting overrrides ground and excited splitting if present.

Here we can see these in action for the D1 states:

atom = "Cs"
g = rq.ground_state(atom, splitting="fs")
e = rq.D1_excited(atom, splitting="hfs")

Cs_cell = rq.Cell(atom, [g,e])
print(Cs_cell.states)

[(n=6, l=0, j=0.5, m_j=-0.5), (n=6, l=0, j=0.5, m_j=0.5), (n=6, l=1, j=0.5, f=3.0, m_f=-3.0), (n=6, l=1, j=0.5, f=3.0, m_f=-2.0), (n=6, l=1, j=0.5, f=3.0, m_f=-1.0), (n=6, l=1, j=0.5, f=3.0, m_f=0.0), (n=6, l=1, j=0.5, f=3.0, m_f=1.0), (n=6, l=1, j=0.5, f=3.0, m_f=2.0), (n=6, l=1, j=0.5, f=3.0, m_f=3.0), (n=6, l=1, j=0.5, f=4.0, m_f=-4.0), (n=6, l=1, j=0.5, f=4.0, m_f=-3.0), (n=6, l=1, j=0.5, f=4.0, m_f=-2.0), (n=6, l=1, j=0.5, f=4.0, m_f=-1.0), (n=6, l=1, j=0.5, f=4.0, m_f=0.0), (n=6, l=1, j=0.5, f=4.0, m_f=1.0), (n=6, l=1, j=0.5, f=4.0, m_f=2.0), (n=6, l=1, j=0.5, f=4.0, m_f=3.0), (n=6, l=1, j=0.5, f=4.0, m_f=4.0)]

Or, equivalently:

Cs_cell2 = rq.Cell(atom, rq.D1_states(atom, g_splitting="fs", e_splitting="hfs"))
print(Cs_cell2.states)

[(n=6, l=0, j=0.5, m_j=-0.5), (n=6, l=0, j=0.5, m_j=0.5), (n=6, l=1, j=0.5, f=3.0, m_f=-3.0), (n=6, l=1, j=0.5, f=3.0, m_f=-2.0), (n=6, l=1, j=0.5, f=3.0, m_f=-1.0), (n=6, l=1, j=0.5, f=3.0, m_f=0.0), (n=6, l=1, j=0.5, f=3.0, m_f=1.0), (n=6, l=1, j=0.5, f=3.0, m_f=2.0), (n=6, l=1, j=0.5, f=3.0, m_f=3.0), (n=6, l=1, j=0.5, f=4.0, m_f=-4.0), (n=6, l=1, j=0.5, f=4.0, m_f=-3.0), (n=6, l=1, j=0.5, f=4.0, m_f=-2.0), (n=6, l=1, j=0.5, f=4.0, m_f=-1.0), (n=6, l=1, j=0.5, f=4.0, m_f=0.0), (n=6, l=1, j=0.5, f=4.0, m_f=1.0), (n=6, l=1, j=0.5, f=4.0, m_f=2.0), (n=6, l=1, j=0.5, f=4.0, m_f=3.0), (n=6, l=1, j=0.5, f=4.0, m_f=4.0)]

2. Decoherence rates in Cell

As we can see above, natural decays based on the atom used are added automatically. While the basic mechanics of how decoherence is handled in Cell is identical to Sensor, there are some differences in what get added to the couplings graph that need to be accounted for. In this section we discuss the details of decoherent couplings in the Cell class, including what is calculated automatically, what isn’t, and how to make it consistent with calculations and experiments outside of rydiqule.

2.1 Natural transition rates

Obviously, in an actual atomic vapor, population in higher energies naturally decays to lower energies. In rydiqule, this is handled automatically via the gamma_transition keyword added to the graph edge. Let us start by recreating the very simple Cell we created in section 1 on the D1 line of Rubidium-85, and examine the edges of the graph before anything else gets added.

g = A_QState(5, 0, 0.5)
e = A_QState(5, 1, 0.5)

Rb_Cell = rq.Cell("Rb85", [g, e])
print(Rb_Cell.couplings.edges(data=True))

[((n=5, l=1, j=0.5), (n=5, l=0, j=0.5), {'gamma_transition': 36.11450417508357, 'label': '((5, 1, 0.5),(5, 0, 0.5))'})]

Here we can indeed see that, unlike in a basic Sensor, the transition rate (in Mrad/s as always) is added from the excited state to the ground state without any extra work on our part. This will be the case no matter how many states are added to the Cell. Let us again examine this attribute but with the full hyperfine manifolds of the D1 line.

atom = "Cs" #Cesium for fun, only one isotope (133) supported by ARC so the isotope number is omitted

D1_hfs_cell = rq.Cell(atom, rq.D1_states(atom, splitting="hfs"))
for s1, s2, value in D1_hfs_cell.couplings.edges(data="gamma_transition"):
    print(f"({s1},{s2}):{value}")

((6, 1, 0.5, f=3.0, m_f=-3.0),(6, 0, 0.5, f=3.0, m_f=-3.0)):5.37334121162191
((6, 1, 0.5, f=3.0, m_f=-3.0),(6, 0, 0.5, f=3.0, m_f=-2.0)):1.7911137372073023
((6, 1, 0.5, f=3.0, m_f=-3.0),(6, 0, 0.5, f=4.0, m_f=-4.0)):16.717061547268166
((6, 1, 0.5, f=3.0, m_f=-3.0),(6, 0, 0.5, f=4.0, m_f=-3.0)):4.1792653868170415
((6, 1, 0.5, f=3.0, m_f=-3.0),(6, 0, 0.5, f=4.0, m_f=-2.0)):0.5970379124024344
((6, 1, 0.5, f=3.0, m_f=-2.0),(6, 0, 0.5, f=3.0, m_f=-3.0)):1.7911137372073023
((6, 1, 0.5, f=3.0, m_f=-2.0),(6, 0, 0.5, f=3.0, m_f=-2.0)):2.388151649609737
((6, 1, 0.5, f=3.0, m_f=-2.0),(6, 0, 0.5, f=3.0, m_f=-1.0)):2.9851895620121716
((6, 1, 0.5, f=3.0, m_f=-2.0),(6, 0, 0.5, f=4.0, m_f=-3.0)):12.537796160451125
((6, 1, 0.5, f=3.0, m_f=-2.0),(6, 0, 0.5, f=4.0, m_f=-2.0)):7.16445494882921
((6, 1, 0.5, f=3.0, m_f=-2.0),(6, 0, 0.5, f=4.0, m_f=-1.0)):1.7911137372073025
((6, 1, 0.5, f=3.0, m_f=-1.0),(6, 0, 0.5, f=3.0, m_f=-2.0)):2.9851895620121702
((6, 1, 0.5, f=3.0, m_f=-1.0),(6, 0, 0.5, f=3.0, m_f=-1.0)):0.5970379124024342
((6, 1, 0.5, f=3.0, m_f=-1.0),(6, 0, 0.5, f=3.0, m_f=0.0)):3.5822274744146045
((6, 1, 0.5, f=3.0, m_f=-1.0),(6, 0, 0.5, f=4.0, m_f=-2.0)):8.955568686036512
((6, 1, 0.5, f=3.0, m_f=-1.0),(6, 0, 0.5, f=4.0, m_f=-1.0)):8.955568686036512
((6, 1, 0.5, f=3.0, m_f=-1.0),(6, 0, 0.5, f=4.0, m_f=0.0)):3.5822274744146045
((6, 1, 0.5, f=3.0, m_f=0.0),(6, 0, 0.5, f=3.0, m_f=-1.0)):3.5822274744146045
((6, 1, 0.5, f=3.0, m_f=0.0),(6, 0, 0.5, f=3.0, m_f=1.0)):3.5822274744146045
((6, 1, 0.5, f=3.0, m_f=0.0),(6, 0, 0.5, f=4.0, m_f=-1.0)):5.970379124024342
((6, 1, 0.5, f=3.0, m_f=0.0),(6, 0, 0.5, f=4.0, m_f=0.0)):9.55260659843895
((6, 1, 0.5, f=3.0, m_f=0.0),(6, 0, 0.5, f=4.0, m_f=1.0)):5.970379124024342
((6, 1, 0.5, f=3.0, m_f=1.0),(6, 0, 0.5, f=3.0, m_f=0.0)):3.5822274744146045
((6, 1, 0.5, f=3.0, m_f=1.0),(6, 0, 0.5, f=3.0, m_f=1.0)):0.5970379124024342
((6, 1, 0.5, f=3.0, m_f=1.0),(6, 0, 0.5, f=3.0, m_f=2.0)):2.9851895620121702
((6, 1, 0.5, f=3.0, m_f=1.0),(6, 0, 0.5, f=4.0, m_f=0.0)):3.5822274744146045
((6, 1, 0.5, f=3.0, m_f=1.0),(6, 0, 0.5, f=4.0, m_f=1.0)):8.955568686036512
((6, 1, 0.5, f=3.0, m_f=1.0),(6, 0, 0.5, f=4.0, m_f=2.0)):8.955568686036512
((6, 1, 0.5, f=3.0, m_f=2.0),(6, 0, 0.5, f=3.0, m_f=1.0)):2.9851895620121716
((6, 1, 0.5, f=3.0, m_f=2.0),(6, 0, 0.5, f=3.0, m_f=2.0)):2.388151649609737
((6, 1, 0.5, f=3.0, m_f=2.0),(6, 0, 0.5, f=3.0, m_f=3.0)):1.7911137372073023
((6, 1, 0.5, f=3.0, m_f=2.0),(6, 0, 0.5, f=4.0, m_f=1.0)):1.7911137372073025
((6, 1, 0.5, f=3.0, m_f=2.0),(6, 0, 0.5, f=4.0, m_f=2.0)):7.16445494882921
((6, 1, 0.5, f=3.0, m_f=2.0),(6, 0, 0.5, f=4.0, m_f=3.0)):12.537796160451125
((6, 1, 0.5, f=3.0, m_f=3.0),(6, 0, 0.5, f=3.0, m_f=2.0)):1.7911137372073023
((6, 1, 0.5, f=3.0, m_f=3.0),(6, 0, 0.5, f=3.0, m_f=3.0)):5.37334121162191
((6, 1, 0.5, f=3.0, m_f=3.0),(6, 0, 0.5, f=4.0, m_f=2.0)):0.5970379124024344
((6, 1, 0.5, f=3.0, m_f=3.0),(6, 0, 0.5, f=4.0, m_f=3.0)):4.1792653868170415
((6, 1, 0.5, f=3.0, m_f=3.0),(6, 0, 0.5, f=4.0, m_f=4.0)):16.717061547268166
((6, 1, 0.5, f=4.0, m_f=-4.0),(6, 0, 0.5, f=3.0, m_f=-3.0)):16.71706154726816
((6, 1, 0.5, f=4.0, m_f=-4.0),(6, 0, 0.5, f=4.0, m_f=-4.0)):9.552606598438953
((6, 1, 0.5, f=4.0, m_f=-4.0),(6, 0, 0.5, f=4.0, m_f=-3.0)):2.388151649609738
((6, 1, 0.5, f=4.0, m_f=-3.0),(6, 0, 0.5, f=3.0, m_f=-3.0)):4.17926538681704
((6, 1, 0.5, f=4.0, m_f=-3.0),(6, 0, 0.5, f=3.0, m_f=-2.0)):12.537796160451125
((6, 1, 0.5, f=4.0, m_f=-3.0),(6, 0, 0.5, f=4.0, m_f=-4.0)):2.388151649609738
((6, 1, 0.5, f=4.0, m_f=-3.0),(6, 0, 0.5, f=4.0, m_f=-3.0)):5.37334121162191
((6, 1, 0.5, f=4.0, m_f=-3.0),(6, 0, 0.5, f=4.0, m_f=-2.0)):4.179265386817042
((6, 1, 0.5, f=4.0, m_f=-2.0),(6, 0, 0.5, f=3.0, m_f=-3.0)):0.5970379124024344
((6, 1, 0.5, f=4.0, m_f=-2.0),(6, 0, 0.5, f=3.0, m_f=-2.0)):7.16445494882921
((6, 1, 0.5, f=4.0, m_f=-2.0),(6, 0, 0.5, f=3.0, m_f=-1.0)):8.955568686036514
((6, 1, 0.5, f=4.0, m_f=-2.0),(6, 0, 0.5, f=4.0, m_f=-3.0)):4.179265386817042
((6, 1, 0.5, f=4.0, m_f=-2.0),(6, 0, 0.5, f=4.0, m_f=-2.0)):2.3881516496097372
((6, 1, 0.5, f=4.0, m_f=-2.0),(6, 0, 0.5, f=4.0, m_f=-1.0)):5.373341211621909
((6, 1, 0.5, f=4.0, m_f=-1.0),(6, 0, 0.5, f=3.0, m_f=-2.0)):1.7911137372073025
((6, 1, 0.5, f=4.0, m_f=-1.0),(6, 0, 0.5, f=3.0, m_f=-1.0)):8.955568686036514
((6, 1, 0.5, f=4.0, m_f=-1.0),(6, 0, 0.5, f=3.0, m_f=0.0)):5.970379124024343
((6, 1, 0.5, f=4.0, m_f=-1.0),(6, 0, 0.5, f=4.0, m_f=-2.0)):5.373341211621909
((6, 1, 0.5, f=4.0, m_f=-1.0),(6, 0, 0.5, f=4.0, m_f=-1.0)):0.5970379124024343
((6, 1, 0.5, f=4.0, m_f=-1.0),(6, 0, 0.5, f=4.0, m_f=0.0)):5.970379124024346
((6, 1, 0.5, f=4.0, m_f=0.0),(6, 0, 0.5, f=3.0, m_f=-1.0)):3.5822274744146045
((6, 1, 0.5, f=4.0, m_f=0.0),(6, 0, 0.5, f=3.0, m_f=0.0)):9.55260659843895
((6, 1, 0.5, f=4.0, m_f=0.0),(6, 0, 0.5, f=3.0, m_f=1.0)):3.5822274744146045
((6, 1, 0.5, f=4.0, m_f=0.0),(6, 0, 0.5, f=4.0, m_f=-1.0)):5.970379124024342
((6, 1, 0.5, f=4.0, m_f=0.0),(6, 0, 0.5, f=4.0, m_f=1.0)):5.970379124024342
((6, 1, 0.5, f=4.0, m_f=1.0),(6, 0, 0.5, f=3.0, m_f=0.0)):5.970379124024343
((6, 1, 0.5, f=4.0, m_f=1.0),(6, 0, 0.5, f=3.0, m_f=1.0)):8.955568686036512
((6, 1, 0.5, f=4.0, m_f=1.0),(6, 0, 0.5, f=3.0, m_f=2.0)):1.7911137372073025
((6, 1, 0.5, f=4.0, m_f=1.0),(6, 0, 0.5, f=4.0, m_f=0.0)):5.970379124024342
((6, 1, 0.5, f=4.0, m_f=1.0),(6, 0, 0.5, f=4.0, m_f=1.0)):0.5970379124024343
((6, 1, 0.5, f=4.0, m_f=1.0),(6, 0, 0.5, f=4.0, m_f=2.0)):5.373341211621909
((6, 1, 0.5, f=4.0, m_f=2.0),(6, 0, 0.5, f=3.0, m_f=1.0)):8.955568686036514
((6, 1, 0.5, f=4.0, m_f=2.0),(6, 0, 0.5, f=3.0, m_f=2.0)):7.16445494882921
((6, 1, 0.5, f=4.0, m_f=2.0),(6, 0, 0.5, f=3.0, m_f=3.0)):0.5970379124024344
((6, 1, 0.5, f=4.0, m_f=2.0),(6, 0, 0.5, f=4.0, m_f=1.0)):5.373341211621909
((6, 1, 0.5, f=4.0, m_f=2.0),(6, 0, 0.5, f=4.0, m_f=2.0)):2.3881516496097372
((6, 1, 0.5, f=4.0, m_f=2.0),(6, 0, 0.5, f=4.0, m_f=3.0)):4.179265386817042
((6, 1, 0.5, f=4.0, m_f=3.0),(6, 0, 0.5, f=3.0, m_f=2.0)):12.537796160451125
((6, 1, 0.5, f=4.0, m_f=3.0),(6, 0, 0.5, f=3.0, m_f=3.0)):4.17926538681704
((6, 1, 0.5, f=4.0, m_f=3.0),(6, 0, 0.5, f=4.0, m_f=2.0)):4.179265386817042
((6, 1, 0.5, f=4.0, m_f=3.0),(6, 0, 0.5, f=4.0, m_f=3.0)):5.37334121162191
((6, 1, 0.5, f=4.0, m_f=3.0),(6, 0, 0.5, f=4.0, m_f=4.0)):2.388151649609738
((6, 1, 0.5, f=4.0, m_f=4.0),(6, 0, 0.5, f=3.0, m_f=3.0)):16.71706154726816
((6, 1, 0.5, f=4.0, m_f=4.0),(6, 0, 0.5, f=4.0, m_f=3.0)):2.388151649609738
((6, 1, 0.5, f=4.0, m_f=4.0),(6, 0, 0.5, f=4.0, m_f=4.0)):9.552606598438953

That’s a lot of couplings! We can already start to see the benefits of using a Cell for calculations involving real atoms, since in a Sensor these would all need to be added manually. As we mentioned previously, rydiqule uses the ARC Rydberg package for all of its under-the-hood calculations in Cell, and for more details about how these numbers are computed you can refer to their documentation. It is worth noting that the numbers returned will not be the same since rydiqule converts all relevant quanties to Mrad/s for internal consistency.

2.2 Natural state lifetimes

As we saw in example 1.1, the nodes of the couplings graph also contain information about the natural lifetime of each state in the gamma_lifetime node attribute. We will recreate that example here and see that, for a simple D1 line, the natural lifetime matches the transition rate from the exited to ground state. As a simple sanity check, we also see that the gamma_lifetime attribute for the ground state is 0, since population will obviously not decay out of \(5S^{1/2}\)

g = A_QState(5, 0, 0.5)
e = A_QState(5, 1, 0.5)

Rb_Cell = rq.Cell("Rb85", [g, e])

print(Rb_Cell.couplings.nodes(data="gamma_lifetime"))
print(Rb_Cell.couplings.edges(data="gamma_transition"))

[((n=5, l=0, j=0.5), 0.0), ((n=5, l=1, j=0.5), 36.11450417508357)]
[((n=5, l=1, j=0.5), (n=5, l=0, j=0.5), 36.11450417508357)]

This is because there is only an single decay path for an atom in the \(5P^{1/2}\) state: straight back to the \(5S^{1/2}\) state. However, if we now imagine that the the first excited state is one which has multiple decay paths, it will no longer be the case that the natural state lifetime matches the decay rate. To see this, let us add the \(6P^{3/2}\) state as the first excited state:

g = A_QState(5, 0, 0.5)
e1 = A_QState(6, 1, 1.5)

Rb_Cell_6S = rq.Cell("Rb85", [g, e1])
print(Rb_Cell_6S.couplings.nodes(data="gamma_lifetime"))
print(Rb_Cell_6S.couplings.edges(data="gamma_transition"))
print(Rb_Cell_6S.decoherence_matrix())

[((n=5, l=0, j=0.5), 0.0), ((n=6, l=1, j=1.5), 8.464336938490163)]
[((n=6, l=1, j=1.5), (n=5, l=0, j=0.5), 1.9967623792713765)]
[[0.         0.        ]
 [8.46433694 0.        ]]

You may notice here that gamma_transition value on the edge is less than the state lifetime, but the decoherence matrix still accounts for the entire state lifetime over all paths. What gives? Where did the rest of that decay rate come from if the decoherence_matrix() method works exactly the same in Sensor as it does in Cell? To answer that question, let us again inspect the graph edges completely.

print(Rb_Cell_6S.couplings.edges(data=True))

[((n=6, l=1, j=1.5), (n=5, l=0, j=0.5), {'gamma_transition': 1.9967623792713765, 'label': '((6, 1, 1.5),(5, 0, 0.5))', 'gamma_mismatch': 6.467574559218787})]

Aha! rydiqule has created an extra attribute on the edge called gamma_mismatch which pretty much does what it sounds like. It is added by rydiqule if the decay rates out of a particular state do not sum to the state lifetime. Since all terms starting with "gamma" on an edge are summed to calculate that particular term in the decoherence matrix, we maintain a transparent convention that preserves the functionality of Sensor.decoherence_matrix without any opaque magic. In fact, if you dig into the source code of Cell, you will not actually find a Cell.decoherence_matrix, it inherits the function directly from Sensor.

Default behavior

To demonstrate this behavior more completely, lets recreate this Cell with all possible decay paths out of \(6P^{3/2}\) and see that there is then no mismatch to account for, all decays sum to the decay rate associated with the lifetime of the state. In this case, there are 4 allowed decays based on selection rules: \(5S^{1/2}\) (the ground state of the atom), \(6S^{1/2}\), \(4D^{3/2}\), and \(4D^{5/2}\)

atom = "Rb85"
g = rq.ground_state(atom)

e_max = A_QState(6, 1, 1.5)

e1 = A_QState(6,0,0.5)
e2 = A_QState(4, 2, 1.5)
e3 = A_QState(4, 2, 2.5)

Rb_Cell_6S_full_decay = rq.Cell(atom, [g, e1, e2, e3, e_max]) 

for s1, s2, gamma_mismatch in Rb_Cell_6S_full_decay.couplings.edges(data="gamma_mismatch"):
    print(f"({s1},{s2}: {gamma_mismatch})")

((6, 0, 0.5),(5, 0, 0.5): 21.750359563468862)
((4, 2, 1.5),(5, 0, 0.5): 11.481506041906623)
((4, 2, 2.5),(5, 0, 0.5): 10.675163106446323)
((6, 1, 1.5),(5, 0, 0.5): None)
((6, 1, 1.5),(6, 0, 0.5): None)
((6, 1, 1.5),(4, 2, 1.5): None)
((6, 1, 1.5),(4, 2, 2.5): None)

We have gamma_mismatch on new edges (those out of states that ought to decay to \(5P\) but don’t because it is not in the system), but importatntly for this demonstration, there are no gamma_mismatch values associated with \(6P^{3/2}\), because the total decay out of that state equals its gamma_lifetime value. We can check this to be sure in the following way:

gamma_lifetime = Rb_Cell_6S_full_decay.couplings.nodes[e_max]["gamma_lifetime"]
print(f"gamma_lifetime: {gamma_lifetime}")

total_decay = 0
for state in [g, e1, e2, e3]:
    total_decay += Rb_Cell_6S_full_decay.couplings.edges[e_max, state]["gamma_transition"]
print(f"total gamma_transition: {total_decay}")

gamma_lifetime: 8.464336938490163
total gamma_transition: 8.464336938490163

You may notice in the above examples that gamma_mismatch is always sent to ground automatically. It is often suffient to assume this since in many systems population will eventually end up back in the ground state regardless of what path it took. Still this is not always the case, so how de we control this behavior? rydiqule has a couple of options that we control using the gamma_mismatch option in the constructor:

1. "ground", the option demonstrated above, is the default. It adds a decoherent coupling from each state to \(nS^{1/2}\), where n is the principle quantum number of the atom specified in the cell. If more that one sublevel exists for this state in the Cell, the decoherence will be divided equally amongst all sublevels. This coupling will be added regardless of whether the transition is dipole-allowed, since it is based on the assumption that all population will eventually decay to ground through some pathway. This option should be avoided for systems containing dark states.

2. "all" is an option that divides the mismatching decay values between all other calculated decay paths. The fraction of gamma_mismatch added to each edge is weighted by the existing calculated natural transition rates. So if state \(|3\rangle\) has dipole-allowed decays to states \(|2\rangle\) and \(|1\rangle\) of 6 and 4 respectively, and a total gamma_lifetime in state \(|3\rangle\) of 15, a gamma_mismatch of 3 and 2 will be added to the (3,2) and (3,1) edges respectively. Note that there must be a dipole-allowed decay for every state in the Cell other than the ground state or else selecting this option will error.

3. "none" does exactly what it sounds like: it will not add a single decay to the system for the computed state lifetime and decay rate mismatch. In this case, it is typically assumed that you will add a decay manually, but this leaves acounting for decay up to you.

In the future, rydiqule might support other options for which there is a compelling use case, but for now, anything besides the ground and all as described above will require manual specification of decays using the none option.

3.3 Transit broadening

As one final note on decoherent transitions, Cell treats transit broadening exactly the same is Sensor, although it does compute a transit broadening rate automatically based on atomic temperature, beam area, and atomic mass (beam is presumbed to be gaussian). So while the value of transit_broadening is computed automatically and stored as an attribute, you need to properly define the ensemble temperature and the optical beam area.

3. Couplings in Cell

Just like most things in Cell, couplings work much the same way as they do in Sensor but with some additional functionality.

3.1 Automatic Quantities

To see the basics of what is added automatically, let us create a simple 2-state system, add a single coupling between them, then inspect the resutling graph.

atom = "Rb85"
[g, e] = rq.D1_states(atom)

Rb_Cell_basic = rq.Cell(atom, [g, e])
Rb_Cell_basic.add_coupling((g, e), rabi_frequency=1, detuning=1, label="laser")

print(Rb_Cell_basic.couplings.edges(data=True))

[((n=5, l=0, j=0.5), (n=5, l=1, j=0.5), {'rabi_frequency': 1, 'detuning': 1, 'phase': 0, 'kvec': (0, 0, 0), 'label': 'laser', 'coherent_cc': 1, 'dipole_moment': 1.7277475900721146, 'q': 0}), ((n=5, l=1, j=0.5), (n=5, l=0, j=0.5), {'gamma_transition': 36.11450417508357, 'label': '((5, 1, 0.5),(5, 0, 0.5))'})]

Here we can see a coupling things. Firstly, the dipole_moment is calculated automatically. While this is not used directly in this case, it can come in handy for time-dependant couplings, and it can be a useful reference for other calculations you may want to do afterwards. Futhermore, if we add a coupling not in the rotating wave approximation, the transition_frequency will be calculated automatically. It is worth noting, however, that for very large transition frequencies, it can take a very long time (often prohibitively long) to solve in the time domain. rydiqule will warn you if you add a coupling without the RWA if the transition frequency is very high. If you are aware of this fact and want to suppress the warning, you can suppressed with warnings.simplefilter('ignore', rq.RWAWarning), but it is often not desired to compute in the time domain when transition frequencies are quite large.

import warnings
warnings.simplefilter("ignore", rq.RWAWarning)

atom = "Rb85"
[g, e] = rq.D1_states(atom)

Rb_Cell_time = rq.Cell(atom, [g, e])
Rb_Cell_time.add_coupling((g, e), rabi_frequency=1, label="laser")

print(Rb_Cell_time.couplings.edges(data=True))

[((n=5, l=0, j=0.5), (n=5, l=1, j=0.5), {'rabi_frequency': 1, 'transition_frequency': 2369435883.8824973, 'phase': 0, 'kvec': (0, 0, 0), 'label': 'laser', 'coherent_cc': 1, 'dipole_moment': 1.7277475900721146, 'q': 0}), ((n=5, l=1, j=0.5), (n=5, l=0, j=0.5), {'gamma_transition': 36.11450417508357, 'label': '((5, 1, 0.5),(5, 0, 0.5))'})]

3.2 Aliasing manifolds with variables

While manifolds can be all defined manually, it is worth a brief sidebar here to discuss how to efficiently do it in Cell when manifolds can get quite large. We can already see one useful trick used in the above cell, by declaring [g, e] = rq.D1_states(atom). In python, lists can be converted to individual single variables something like [a,b,c] = my_list as long as my_list has the number of elements that are upacked in the list defined.

Now we note that helper functions like D1_states and expand_qnums, which we introduce below, return lists of states. This allows for easy aliasing of states to the extent that we would like to apply them in couplings. We already showed this functionality to in the example above with [g, e] = rq.D1_states(atom), but we can show some more involved examples here. Note that just like the constructor of Cell, states are passed to expand_qnums as a list.

While contrived, this example will show the utility of this approach well. Here we add a coupling from the ground state to the excited manifold, but add decoherences only from the \(m_j = \pm0.5\) states back to ground.

atom = "Rb85"
g = rq.ground_state(atom, splitting="fs")
e_fs = rq.D2_excited(atom, splitting="fs") #4 states total
print(f"states: {g, e_fs}\n")
[e1, e2, e3, e4] = rq.expand_qnums([e_fs])
print(f"m_j = +\-0.5 states: {e2, e3}\n")

my_cell = rq.Cell(atom, [g, e_fs])
my_cell.add_coupling((g,e1), rabi_frequency=1, detuning=0, label="laser")
my_cell.add_decoherence((e2, g), 1, label="foo")
my_cell.add_decoherence((e3, g), 1, label="bar")

print("couplings")
for edge in my_cell.couplings.edges(data=True):
    print(edge)

print(f"gamma matrix: \n {my_cell.decoherence_matrix()}")

states: ((n=5, l=0, j=0.5, m_j='all'), (n=5, l=1, j=1.5, m_j='all'))

m_j = +\-0.5 states: ((n=5, l=1, j=1.5, m_j=-0.5), (n=5, l=1, j=1.5, m_j=0.5))

couplings
((n=5, l=1, j=1.5, m_j=-1.5), (n=5, l=0, j=0.5, m_j=-0.5), {'gamma_transition': 38.11316014416465, 'label': '((5, 1, 1.5, m_j=-1.5),(5, 0, 0.5, m_j=-0.5))'})
((n=5, l=1, j=1.5, m_j=-0.5), (n=5, l=0, j=0.5, m_j=-0.5), {'gamma_transition': 25.408773429443098, 'label': '((5, 1, 1.5, m_j=-0.5),(5, 0, 0.5, m_j=-0.5))', 'gamma_foo': 1.0})
((n=5, l=1, j=1.5, m_j=-0.5), (n=5, l=0, j=0.5, m_j=0.5), {'gamma_transition': 12.704386714721549, 'label': '((5, 1, 1.5, m_j=-0.5),(5, 0, 0.5, m_j=0.5))', 'gamma_foo': 1.0})
((n=5, l=1, j=1.5, m_j=0.5), (n=5, l=0, j=0.5, m_j=-0.5), {'gamma_transition': 12.704386714721549, 'label': '((5, 1, 1.5, m_j=0.5),(5, 0, 0.5, m_j=-0.5))', 'gamma_bar': 1.0})
((n=5, l=1, j=1.5, m_j=0.5), (n=5, l=0, j=0.5, m_j=0.5), {'gamma_transition': 25.408773429443098, 'label': '((5, 1, 1.5, m_j=0.5),(5, 0, 0.5, m_j=0.5))', 'gamma_bar': 1.0})
((n=5, l=1, j=1.5, m_j=1.5), (n=5, l=0, j=0.5, m_j=0.5), {'gamma_transition': 38.11316014416465, 'label': '((5, 1, 1.5, m_j=1.5),(5, 0, 0.5, m_j=0.5))'})
gamma matrix: 
 [[ 0.          0.          0.          0.          0.          0.        ]
 [ 0.          0.          0.          0.          0.          0.        ]
 [38.11316014  0.          0.          0.          0.          0.        ]
 [26.40877343 13.70438671  0.          0.          0.          0.        ]
 [13.70438671 26.40877343  0.          0.          0.          0.        ]
 [ 0.         38.11316014  0.          0.          0.          0.        ]]

Breaking states out like this is not always necessary, but it can be useful in certain cases where you would like to only couple certain sublevels. Rather that retyping out A_QState(n,l,j,...) every time you reference it, you can assign it to a single terse variable without actually typing it out manually at all.

Often you only care about couplings between entire manifolds anyway, but this is a useful trick to have in your pocket.

3.3 Alternate specifications of rabi_frequency

Since Cell tries to support an interface to atomic physics that is a more close to real-life experiments, it has a couple of ways to speciy the power of a field beyond the rabi_frequency that is used in Sensor. As we have already demonstrated, you can define the rabi frequency of a coupling just fine (a Cell is just a Sensor after all), but Cell does also introduce a couple more options. All totalled, the three options for rabi_frequency specification are as follows:

1. Base specification of rabi_frequency just as in Sensor.

2. Specification of electric field, in V/m via the the e_field arguments. rydiqule will calculate the rabi_frequency based on the computed dipole moment.

3. Specification of both beam power, in watts and \(\frac{1}{e^2}\) beam waist (radius) in meters via the beam_power and beam_waist arguments. In all cases, the arguments are provided as optional keyword arguments to the add_coupling function. These options are mutually exclusive. Also, none of this information will be stored on the graph directly - it will be converted to rabi_frequency, which is the only such quantity that is stored on the graph. We briefly demonstrate all of this below by defining a simple 3-level vee scheme.

atom = "Rb85"
g = rq.ground_state(atom)
D1_e = rq.D1_excited(atom)
D2_e = rq.D2_excited(atom)

RbCell_multi_rabi = rq.Cell(atom, [g, D1_e, D2_e])
RbCell_multi_rabi.add_coupling((g, D1_e), detuning=1, e_field=1)
RbCell_multi_rabi.add_coupling((g, D2_e), detuning=1, beam_power=1, beam_waist=0.01)

#rabi_frequency is on all edges despite not being directly specified
print(RbCell_multi_rabi.couplings.edges.data("rabi_frequency"))
print(RbCell_multi_rabi.couplings.edges.data("e_field"))
print(RbCell_multi_rabi.couplings.edges.data("beam_power"))

[((n=5, l=0, j=0.5), (n=5, l=1, j=0.5), 0.13890429081447606), ((n=5, l=0, j=0.5), (n=5, l=1, j=1.5), 429.7419875788027), ((n=5, l=1, j=0.5), (n=5, l=0, j=0.5), None), ((n=5, l=1, j=1.5), (n=5, l=0, j=0.5), None)]
[((n=5, l=0, j=0.5), (n=5, l=1, j=0.5), None), ((n=5, l=0, j=0.5), (n=5, l=1, j=1.5), None), ((n=5, l=1, j=0.5), (n=5, l=0, j=0.5), None), ((n=5, l=1, j=1.5), (n=5, l=0, j=0.5), None)]
[((n=5, l=0, j=0.5), (n=5, l=1, j=0.5), None), ((n=5, l=0, j=0.5), (n=5, l=1, j=1.5), None), ((n=5, l=1, j=0.5), (n=5, l=0, j=0.5), None), ((n=5, l=1, j=1.5), (n=5, l=0, j=0.5), None)]

3.4 Clebsch-Gordon coefficients

When defining couplings over manifolds, it is important to consider the Clebsch-Gordon coefficients (which we will refer to as CGC), which define the relating weighting of dipole moments between different sublevels. Recall that in Sensor, this is handled using the coherent_cc optional keyword argument in add_coupling. In Cell, the coherent_cc is filled in automatically using ARC’s functions for the spherical component of the dipole matrix element. This section is really only intended to demonstrate that coherent_cc is calculated automatically in cell. For a more detailed description of how rydiqule handles coupling between manifolds of states, check the physics documentation. Here we just create a simple Cell on the D1 transition with fine structure splitting, and show that coherent_cc is populated on the graph even if we don’t specify it explicitly.

atom = "Rb85"
g = rq.ground_state(atom, splitting="fs")
e = rq.D1_excited(atom, splitting="fs")
[g1, g2] = rq.expand_qnums([g])
[e1, e2] = rq.expand_qnums([e])

RbCell_cc = rq.Cell(atom, [g,e])
RbCell_cc.add_coupling((g,e), rabi_frequency=1, detuning=1, label="D1")
print(RbCell_cc.couplings.edges.data("coherent_cc"))

[((n=5, l=0, j=0.5, m_j=-0.5), (n=5, l=1, j=0.5, m_j=-0.5), -0.816496580927726), ((n=5, l=0, j=0.5, m_j=0.5), (n=5, l=1, j=0.5, m_j=0.5), 0.816496580927726), ((n=5, l=1, j=0.5, m_j=-0.5), (n=5, l=0, j=0.5, m_j=-0.5), None), ((n=5, l=1, j=0.5, m_j=-0.5), (n=5, l=0, j=0.5, m_j=0.5), None), ((n=5, l=1, j=0.5, m_j=0.5), (n=5, l=0, j=0.5, m_j=-0.5), None), ((n=5, l=1, j=0.5, m_j=0.5), (n=5, l=0, j=0.5, m_j=0.5), None)]

It is worth noting that rydiqule is making a very specific choice here regarding its convention regarding coupling coefficients. We choose to use the spherical dipole moment, since it is at least proportial to the Clebsch-Gordon coefficients and should ultimately produce the correct rabi frequency in the Hamiltonian. For consistency, we do not allow for manual specification of coupling coefficients in Cell, so make sure you are familiarized with the documentation above to ensure that the numbers in your simulation are correct.

3.5 Doppler Shifts

Recall that when we want to define a system with doppler broadening in a Sensor, we define the vP parameter either in the constructor or by accessing the attribute directly after construction, and then specifying the kvec parameter in couplings. In a Cell, however, the magnitude of the kvec parameter is fixed by the transition frequency, and the most probable speed is already set by the temperature. So when adding a coupling with doppler broadening in a Cell, we set only the newly-introduced kunit when couplings are defined. This is a tuple representing an (x, y, z) unit vector representing the direction of the plane wave of the field. Other than that, solving with doppler broadening in a Cell is identical to doing so in Sensor, and we can otherwise do everything the same way.

Note that for consistency with Sensor, the same kvec attribute is added to the graph edge, it is just computed automatically from kunit.

atom = "Rb85"
g = rq.ground_state(atom, splitting="fs")
e = rq.D1_excited(atom, splitting="fs")
[g1, g2] = rq.expand_qnums([g])
[e1, e2] = rq.expand_qnums([e])

RbCell_dop = rq.Cell(atom, [g,e])
RbCell_dop.add_coupling((g,e), rabi_frequency=1, detuning=1, label="D1", kunit=(1,0,0))
print(RbCell_dop)

print(rq.solve_steady_state(RbCell_dop, doppler=True).rho)

<class 'rydiqule.cell.Cell'> object with 4 states and 2 coherent couplings.
States: [(n=5, l=0, j=0.5, m_j=-0.5), (n=5, l=0, j=0.5, m_j=0.5), (n=5, l=1, j=0.5, m_j=-0.5), (n=5, l=1, j=0.5, m_j=0.5)]
Coherent Couplings: 
    ((5, 0, 0.5, m_j=-0.5),(5, 1, 0.5, m_j=-0.5)): {rabi_frequency: 1, detuning: 1, phase: 0, kvec: <parameter with 3 values>, label: D1_0, coherent_cc: -0.816496580927726, dipole_moment: -1.7277475900721146, q: 0}
    ((5, 0, 0.5, m_j=0.5),(5, 1, 0.5, m_j=0.5)): {rabi_frequency: 1, detuning: 1, phase: 0, kvec: <parameter with 3 values>, label: D1_1, coherent_cc: 0.816496580927726, dipole_moment: 1.7277475900721146, q: 0}
Decoherent Couplings:
    ((5, 1, 0.5, m_j=-0.5),(5, 0, 0.5, m_j=-0.5)): {gamma_transition: 12.03816805836119}
    ((5, 1, 0.5, m_j=-0.5),(5, 0, 0.5, m_j=0.5)): {gamma_transition: 24.07633611672238}
    ((5, 1, 0.5, m_j=0.5),(5, 0, 0.5, m_j=-0.5)): {gamma_transition: 24.07633611672238}
    ((5, 1, 0.5, m_j=0.5),(5, 0, 0.5, m_j=0.5)): {gamma_transition: 12.03816805836119}
Energy Shifts:
    None
[ 0.00000000e+00 -1.09387364e-07  0.00000000e+00  0.00000000e+00
  4.97771560e-01  0.00000000e+00  1.09387364e-07 -1.86770118e-04
  0.00000000e+00  4.22260158e-06  0.00000000e+00  0.00000000e+00
  1.86770118e-04  0.00000000e+00  4.22260158e-06]

4. Solving systems defined in Cell

Fortunately, this will be a straightforward section. Since Cell inherits Sensor, the mechanics of solving are identical between the two classes. There are a couple of additional considerations about automatically calculated quantities that we will go over in this section.

4.1 Just like Sensor!

Just like in Sensor, calling rq.solve_steady_state or rq.solve_time will produce a solution object containing information pertinent to the solve, as we demonstrate below.

atom = "Rb85"
g = rq.ground_state(atom)
e1 = rq.D2_excited(atom)
e2 = A_QState(6, 2, 2.5)

det = np.linspace(-1,1,11)
rabi = np.linspace(0.1, 1, 10)

RbCell_time = rq.Cell(atom, [g, e1, e2])
RbCell_time.add_coupling((g,e1), rabi_frequency=1, detuning=det, label="D1")
RbCell_time.add_coupling((e1, e2), rabi_frequency=rabi, detuning=1, label="upper")

sol = rq.solve_steady_state(RbCell_time)
print(sol.rho.shape)
print(sol.axis_labels)

(11, 10, 8)
['D1_detuning', 'upper_rabi_frequency', 'density_matrix']

4.2 Observables in Cell

There are a couple other additions of note in Cell that are not present in Sensor. Specifically, some automatic calculations that, in Sensor, would be specified manually. Recall that the kappa and eta quantities, when relevant for calculating observables, would need to be spcecified as defined in the API docs. Since all the relevant parts are already in Cell, it cannot be specified in Cell, it must be computed automatically.

[g,e] = rq.D1_states(5)
RbCell_kappa_eta = rq.Cell("Rb85", [g, e])
RbCell_kappa_eta.add_coupling((g,e), rabi_frequency=1, detuning=0)

print(f"kappa: {RbCell_kappa_eta.kappa}")
print(f"eta: {RbCell_kappa_eta.eta}")

kappa: 14251202077.430082
eta: 0.9529640302753843