from pylgs.lgssystem import LGSSystem
from pylgs.velocitygroups import VelocityGroups
from pymor.algorithms.simplify import contract
from pymor.algorithms.to_matrix import to_matrix
import numpy as npTechnical background
The Liouville-von Neumann equation
The state of an atomic ensemble is described by the density matrix \(\rho_{nm}\), where the indices \(n\) and \(m\) represent attainable eigenstates of the atomic system. The Hamiltonian evolution of the density matrix is given by the Liouville–von Neumann equation \[i\hbar\dot\rho = [H,\rho] \tag{1}\] with additional phenomenological terms added to the right-hand side to account for relaxation and repopulation processes resulting from external interactions.
Writing the \(N_\rho\) matrix elements of \(\rho\) as a vector \(\rho_i\) with a single index \(i\), the density-matrix evolution equation can be written as a linear matrix equation \[\dot\rho_i = A_{ij}\rho_j + b_i, \tag{2}\] where \(\rho\), \(A\), and \(b\) may all depend on time. (Here and below we use Einstein summation notation, in which repeated indices are assumed to be summed over.)
The laser guide star system
Velocity dependence and discretization
Because of the Doppler shift for atoms with nonzero velocity along the direction of the laser beam, it is important to distinguish between atoms with different longitudinal velocities – \(\rho(v)\), \(A(v)\), and \(b(v)\) are all functions of longitudinal velocity \(v\). The constant term \(b(v)\) describes the effect of transit repopulation, while the \(A(v)\) includes the effect of several physical processes, of which Doppler shift, velocity-changing collisions (vcc), and recoil depend on velocity, each in a different manner.
The Doppler shift of the pump light is proportional to \(v\), so the term describing the Doppler shift can be written \[\dot\rho_i(v)|_\text{dop}=vA^\text{dop}_{ij}\rho_j(v).\] Velocity-changing collisions are assumed to change the velocity of an atom to any other velocity with probability given by the Maxwell-Boltzmann distribution \(n(v)\), so this term is proportional to \(n(v)\) and is multiplied by the integral of \(\rho\) over all velocities, \(\int\rho_j(v')dv'\): \[\dot\rho_i(v)|_\text{vcc}=A^\text{vcc}_{ij}n(v)\int\rho_j(v')dv'.\] Atoms undergoing recoil have their velocities increased, transferring them to a neighboring velocity group, so this term is proportional to the derivative of \(\rho(v)\): \[\dot\rho_i(v)|_\text{rec}=A^\text{rec}_{ij}\frac{d\rho_j}{dv}.\] Atoms entering the illuminated region have the Maxwell-Boltzmann velocity distribution, so this term is proportional to \(n(v)\): \[\dot\rho_i(v)|_\text{trans}=n(v)b_i.\] Finally, there is also a term that is independent of the atomic velocity: \[\dot\rho_i(v)|_\text{ind}=A^\text{ind}_{ij}\rho_j(v).\]
Corresponding equations for the density matrices \(\rho_{vi}\) for a discretized velocity distribution with \(N_v\) velocity bins can be obtained by integrating the terms over each bin with mean velocities \(v\), widths \(\Delta v\), and fractional densities \(n_v\Delta v\). This results in the terms \[\dot\rho_{vi}|_\text{dop}=vA^\text{dop}_{ij}\rho_{vj}=(v\delta_{vv'})A^\text{dop}_{ij}\rho_{v'j}, \] \[\dot\rho_{vi}|_\text{vcc}=A^\text{vcc}_{ij}n_v\sum_{v'}\rho_{v'j}=(n_v1_{v'})A^\text{vcc}_{ij}\rho_{v'j},\] \[\dot\rho_{vi}|_\text{rec}=A^\text{rec}_{ij}\frac{\rho_{v+1,j}-\rho_{vj}}{\Delta v}=\frac{\delta_{v+1,v'}-\delta_{vv'}}{\Delta v}A^\text{rec}_{ij}\rho_{v'j},\] \[\dot\rho_{vi}|_\text{trans}=n_vb_i,\] \[\dot\rho_{vi}|_\text{ind}=A^\text{ind}_{ij}\rho_{vj}=\delta_{vv'}A^\text{ind}_{ij}\rho_{v'j},\] where we have inserted operators on velocity space (\(\delta_{vv'}\) is the Kronecker delta function, i.e., identity matrix, and \(1_{v'}\) is 1 for all values of \(v'\)) to make it explicit that we now have a linear equation for the \((N_v\times N_\rho)\) matrix \(\rho_{vi}\) in terms of four-dimensional \((N_v\times N_v \times N_\rho\times N_\rho)\) matrices. The combined equation now reads \[ \dot\rho_{vi} = \left[ \delta_{vv'}A^\text{ind}_{ij} + (v\delta_{vv'})A^\text{dop}_{ij} + (n_v1_{v'})A^\text{vcc}_{ij} + \frac{\delta_{v+1,v'}-\delta_{vv'}}{\Delta v}A^\text{rec}_{ij} \right]\rho_{v'j} + n_vb_i. \tag{3}\]
Representation as a bilinear system
To perform model order reduction on the density-matrix degrees of freedom of the system, we write the discretized evolution equation for one particular velocity group as a bilinear system of the form
\[ \begin{split} E x'(t) & = A x(t) + \sum_{i = 1}^m{N_i x(t) u_i(t)} + B u(t), \\ y(t) & = C x(t) + D u(t), \end{split} \]
where \(u(t)=(u_1(t),\ldots,u_m(t))\) is a vector of inputs.
Splitting the terms into those that depend on atoms of velocity \(v\) and those that depend on all other atoms of velocity \(v'\), we have
\[ \dot\rho_{vi} = \left[ A^\text{ind}_{ij} + vA^\text{dop}_{ij} + n_vA^\text{vcc}_{ij} - \frac{1}{\Delta v}A^\text{rec}_{ij} \right]\rho_{vj} +\sum_{v'\ne v} \left[ n_vA^\text{vcc}_{ij} + \frac{\delta_{v+1,v'}}{\Delta v}A^\text{rec}_{ij} \right]\rho_{v'j} + n_vb_i. \]
\[ \dot\rho_{vi} = \left[ A^\text{ind}_{ij} + n_vA^\text{vcc}_{ij} + vA^\text{dop}_{ij} - \frac{1}{\Delta v}A^\text{rec}_{ij} \right]\rho_{vj} + n_vb_i +A^\text{vcc}_{ij}\sum_{v'\ne v} n_v\rho_{v'j} + \frac{1}{\Delta v}A^\text{rec}_{ij}\rho_{v+1,j}. \]
Discard the vcc terms and assume that \(A^\text{rec}_{ij}\rho_{v+1,j}\) is proportional to \(b_i\)
\[x'= A^\text{ind}x + \sum_i N_ixu_i + Bu. \]
with \(x=\rho_{v}\) and \[u=(n_v + \frac{1}{\Delta v} \sum_{excited states}\rho_{v+1,j}, v,-\frac{1}{\Delta v}),\] \[A=A^\text{ind},\] \[N=(0, A^\text{vcc}, A^\text{dop}, A^\text{rec}),\] \[B=(b, 0, 0)\]
Representation in the pyMOR framework
In the pyMOR framework, vectors of independent variables such as \(\rho_{vi}\) are represented by VectorArray objects, and the matrices that connect them are represented by Operator objects. The pyLGS LGSSystem object can generate the objects corresponding to the terms in Equation 3.
Import the package:
Set numpy print options:
np.set_printoptions(formatter={'float': lambda x: f'{x:^ 8.2}' if x else f'{0:^ 8}'}, linewidth=140)For simplicity, select a toy LGS atomic system with no angular momentum and fix all parameters:
lgs = LGSSystem(
'NaD1_Toy',
fixed_params={'IntensitySI1': 46.,
'DetuningHz1': -6.268e8,
'LaserWidthHz1': 10.0e6,
'BeamTransitRatePerS': 131.944,
'VccRatePerS': 28571.42,
'TemperatureK': 185.0,
'RecoilParameter': 1,
}
)The \(A^\text{ind}_{ij}\) operator is given by lgs.A_ind:
lgs.A_ind
A_ind{DM(4)} → {DM(4)}
The toy atomic system has only two states and four density-matrix elements, so \(A^\text{ind}_{ij}\) is a \(4\times4\) matrix. We can see it explicitly using to_numpy():
lgs.A_ind.matrix<xarray.DataArray 'A_ind' (Density matrix (range): 4, Density matrix (source): 4)> Size: 264B <COO: shape=(4, 4), dtype=float64, nnz=11, fill_value=0.0> Coordinates: * Density matrix (range) (Density matrix (range)) <U50 800B 'ρ<sub>Re, 3S... * Density matrix (source) (Density matrix (source)) <U50 800B 'ρ<sub>Re, 3...
lgs.A_ind.to_numpy()array([[-2.9e+04, 0 , 3.7e+07, 6.1e+07],
[ 0 , -6.2e+07, -3.9e+09, 0 ],
[-1.9e+07, 3.9e+09, -6.2e+07, 1.9e+07],
[ 0 , 0 , -3.7e+07, -6.1e+07]])We can similarly view \(A^\text{dop}_{ij}\):
lgs.A_dop.to_numpy()array([[ 0 , 0 , 0 , 0 ],
[ 0 , 0 , -3.9e+09, 0 ],
[ 0 , 3.9e+09, 0 , 0 ],
[ 0 , 0 , 0 , 0 ]])\(A^\text{vcc}_{ij}\):
lgs.A_vcc.to_numpy()array([[ 2.9e+04, 0 , 0 , 0 ],
[ 0 , 0 , 0 , 0 ],
[ 0 , 0 , 0 , 0 ],
[ 0 , 0 , 0 , 2.9e+04]])\(A^\text{rec}_{ij}\):
lgs.A_rec.to_numpy()array([[ 0 , 0 , 0 , -4.9e+03],
[ 0 , 0 , 0 , 0 ],
[ 0 , 0 , 0 , 0 ],
[ 0 , 0 , 0 , 0 ]])And \(b_{i}\):
lgs.b.to_numpy()array([[ 1.3e+02],
[ 0 ],
[ 0 ],
[ 0 ]])To form the velocity-space operators from Equation 3, we first define a set of velocity groups, here three evenly spaced groups:
vg = VelocityGroups(3)The explicit form of \(\delta_{vv'}\) is
vg.identity().to_numpy()array([[ 1.0 , 0 , 0 ],
[ 0 , 1.0 , 0 ],
[ 0 , 0 , 1.0 ]])\(v\delta_{vv'}\):
vg.velocity_diagonal().to_numpy()array([[ -2.0 , 0 , 0 ],
[ 0 , 0 , 0 ],
[ 0 , 0 , 2.0 ]])\(n_v1_{v'}\):
vg.n_times_1().to_numpy()array([[ 0.079 , 0.079 , 0.079 ],
[ 0.84 , 0.84 , 0.84 ],
[ 0.079 , 0.079 , 0.079 ]])And \((\delta_{v+1,v'}-\delta_{vv'})/\Delta v\):
vg.drho_dv().to_numpy()array([[ 0.5 , 0 , 0 ],
[ -0.5 , 0.5 , 0 ],
[ 0 , -0.5 , 0.5 ]])The combined operator \(A\), including both density-matrix and velocity-group factors, can be constructed with
A = contract(lgs.operator(vg))
A
X{AV(3) ⨉ DM(4)} → {AV(3) ⨉ DM(4)}
Internally, it is represented by a 4-dimensional xarray DataArray:
A.matrix<xarray.DataArray (Atomic velocity (range): 3, Atomic velocity: 3,
Density matrix (range): 4, Density matrix (source): 4)> Size: 4kB
<COO: shape=(3, 3, 4, 4), dtype=float64, nnz=102, fill_value=-0.0>
Coordinates:
* Atomic velocity (range) (Atomic velocity (range)) float64 24B -2.0 0.0 2.0
* Atomic velocity (Atomic velocity) float64 24B -2.0 0.0 2.0
* Density matrix (range) (Density matrix (range)) <U50 800B 'ρ<sub>Re, 3S...
* Density matrix (source) (Density matrix (source)) <U50 800B 'ρ<sub>Re, 3...We can also express it as a \(12\times12\) block matrix with the \(4\times4\) elements of each block referring to density-matrix elements, and each of the \(3\times3\) array of blocks referring to a pair of velocity groups:
to_matrix(A).toarray()array([[ 2.6e+04, 0 , -3.7e+07, -6.1e+07, -2.2e+03, 0 , 0 , 0 , -2.2e+03, 0 , 0 , 0 ],
[ 0 , 6.2e+07, -3.9e+09, 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ],
[ 1.9e+07, 3.9e+09, 6.2e+07, -1.9e+07, 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ],
[ 0 , 0 , 3.7e+07, 6.1e+07, 0 , 0 , 0 , -2.2e+03, 0 , 0 , 0 , -2.2e+03],
[-2.4e+04, 0 , 0 , -2.5e+03, 4.6e+03, 0 , -3.7e+07, -6.1e+07, -2.4e+04, 0 , 0 , 0 ],
[ 0 , 0 , 0 , 0 , 0 , 6.2e+07, 3.9e+09, 0 , 0 , 0 , 0 , 0 ],
[ 0 , 0 , 0 , 0 , 1.9e+07, -3.9e+09, 6.2e+07, -1.9e+07, 0 , 0 , 0 , 0 ],
[ 0 , 0 , 0 , -2.4e+04, 0 , 0 , 3.7e+07, 6.1e+07, 0 , 0 , 0 , -2.4e+04],
[-2.2e+03, 0 , 0 , 0 , -2.2e+03, 0 , 0 , -2.5e+03, 2.6e+04, 0 , -3.7e+07, -6.1e+07],
[ 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 6.2e+07, 1.2e+10, 0 ],
[ 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 1.9e+07, -1.2e+10, 6.2e+07, -1.9e+07],
[ 0 , 0 , 0 , -2.2e+03, 0 , 0 , 0 , -2.2e+03, 0 , 0 , 3.7e+07, 6.1e+07]])The expanded form is used when the solver needs an explicit form of the matrix for sparse ILU factorization.