Tutorial

Introduction

In this page we discuss the basic approach implemented in pyRET and discuss a few examples of how to use the various packages.

Let us say that we have to compute the photon mediated near field energy transfer between a source (labeled as S) and the absorber (labeled as A). In the dilute limit, the isolated emitters comprises discrete energy levels with corresponding energy eigenstates- ground state \(\mid \Phi_{GS}^{(S/A)} \rangle\) and the excited state \(\mid \Phi_{ES}^{(S/A)} \rangle\). These states are assumed to be multireference in general, and thus can be expressed as a linear combination of single Slater determinants. For example,

(1)\[\begin{split}\langle \mid \Phi_{GS}^{(S)} \rangle = \sum_{\substack{j \in \mathrm{occ} \\ i \in \mathrm{unocc}}} \alpha_{ij}^{(GS)} c_j^{(S/A)\dagger} c_i^{(S)} \langle \mid D^{(S)} \rangle\end{split}\]

and similar expressions would exist for the other states. Here \(\langle \mid D^{(S)} \rangle\) is a single Slater determinant built from filling out the first N orbitals in the single electron picture. For details, please see Appendix C and D of Ref. [1].

Further, we aim to compute the matrix element corresponding to transitions between the multireference states with the interaction Hamiltonian \(H_{int}\) [1].

(2)\[H_{int} = \sum_{E = S, A} \sum_{i=1}^{N_E}\left[-\frac{e\hat{\textbf{p}}_i\cdot\hat{\textbf{A}}}{2m_0}-\frac{e\hat{\textbf{A}}\cdot\hat{\textbf{p}}_i}{2m_0}+e\hat{A}_0 +g\frac{e\hbar}{2m_0}\hat{\bf{\sigma}}_i\cdot\bar\nabla\times\hat{\textbf{A}}\right].\]

We employ the Slater-Condon rule to reduce the matrix element between Slater determinant to matrix element between orbitals. For example,

(3)\[\begin{split}\langle \Phi_{ES}^{(S)} \mid H_{\mathrm{int}} \mid \Phi_{GS}^{(S)}, 1_{k,\alpha s} \rangle =\sum_{\substack{j \in \mathrm{occ} \\ i \in \mathrm{unocc}}} \alpha^{*}_{ij} \langle \phi^{(S)}_{j} \mid H_{\mathrm{int}} \mid \phi^{(S)}_{i}, 1_{k,\alpha s} \rangle\end{split}\]

The single particle orbitals here are \(\left|\phi^{(S)}_{i}\right\rangle\). This code provides an example of how we extract and store these single particle wavefunctions calculated using the Quantum Espresso Density Functional Theory code.

The matrix element corresponding to the photon emission process to a photon mode of energy \(\hbar c k\) and \(\alpha\) denoting all other degrees of freedom (e.g. in this case orbital angular momentum L, total z-projected angular momentum M, and parity P) can be expressed as,

(4)\[v_{k,\alpha}^{(S)} \sqrt{\Delta k} = \langle \Phi_{GS}^{(S)}, 1_{k,\alpha} \mid H_{\mathrm{int}} \mid \Phi_{ES}^{(S)} \rangle\]

and the matrix element corresponding to the photon absorption process can be expressed as,

(5)\[v_{k,\alpha}^{(A)} \sqrt{\Delta k} = \langle \Phi_{ES}^{(A)} \mid H_{\mathrm{int}} \mid \Phi_{GS}^{(A)}, 1_{k,\alpha} \rangle.\]

Here \(\Delta k\) is the mode width parameter used to normalize the photon modes. For details, see Ref. [1].

Once the matrix elements are computed, we then computed the overall energy transfer process matrix element as,

(6)\[M = \frac{i\pi n_i}{\hbar c}\sum_{k,\alpha} v_{k,\alpha}^{(S)} v_{k,\alpha}^{(A)}.\]

Once M is computed, the energy transfer can be computed under second order perturbation using Fermi’s golden rule as outlined in Ref. [1].

With this background, now to start the proess we need the orbital wavefuntions that constitute the multireference states. In the following section, we discuss how the wavefunctions are handled in pyRET.

Loading Wavefunctions

The wavefunction object in pyRET is designed to contain three categories of information:

(1) Plane wave expansion

A plane wave expansion of the wavefunction, as typically done in DFT codes such as Quantum ESPRESSO. In plane wave basis, the single electron orbitals are defined as

(7)\[\phi(r) = \frac{1}{\sqrt{V}} \sum_{\mathbf{G}} c^{-}_{n\mathbf{k}}(\mathbf{G})e^{i(\mathbf{k}+\mathbf{G}) \cdot \mathbf{r}}\]

The plane wave vectors \(\mathbf{G}\) and the expansion coefficients \(c\) are stored as .dat or .hdf5 files outputed by Quantum Espresso. pyRET provides modules that can read these plane wave basis and store the coefficients in an object of the WFunction class. The plane wave expansion is used to directly compute the k-operated wavefunction values, which are used in the computation of the matrix elements of the interaction Hamiltonian. For example,

(8)\[\begin{split}\begin{aligned} k_x \phi(r) &= -i \frac{\partial}{\partial x} \phi(r) \\ &= \frac{1}{\sqrt{V}} \sum_{\mathbf{G}} c^{-}_{n\mathbf{k}}(\mathbf{G}) (k_x+G_x) e^{i(\mathbf{k}+\mathbf{G}) \cdot \mathbf{r}} \end{aligned}\end{split}\]

(2) Spherical harmonic expansion

Alternatively, pyRET can operate the WFunction class in the form of a multipolar expansion of the wavefunction using the Scalar_field class where the angular components are expressed using the MP_SA (multipolar-scalar-angular) class and the radial function can be chosen to be either spherical Bessel functions or atomic orbital like radial dependence. The last feature helps in creating atomic orbital basis within pyRET.

(3) Real space grid

A real space grid of \(\mid \phi(r)\rangle\), \(k_x \mid \phi(r)\rangle\), \(k_y \mid \phi(r)\rangle\), and \(k_z \mid \phi(r)\rangle\). These k-operated values are directly computed from the plane wave expansion, and thus free from any numerical differentiation errors in computing the spatial derivatives. The calculated position grid are stored as an object of the WFunctions_data class defined in the data_classes.py module. For more details please see the examples in the examples directory.

Define emitters and their positions

The next step is to define the emitter types and emitter positions corresponding to the system we want to simulate. pyRET stores the information on the emitters (source and absorber) and their relative positions as objects of classes Emitters_data and Positions_data. The emitters can be defined as dipolar emitters (electric dipole or magnetic dipole), or realistic emitters the orbitals of which are defined from first principles. The objects of the Emitters_data class stores as a dictionary all the defined single-orbital transitions, radiative lifetimes, energies, and other optional properties (like transition dipole moment for ideal dipole-like emitters) where the keys of the dictionary represents the emitter-id (represented as a string). Each emitter in the system of emitters possess an unique id. However, emitters with different emitter id can have same emitter types (also representd as a string). The Positions_data class objects store the relative positions between all possible pairs of emitters. Please refer to the detailed documentation and the example scripts for more details on usage.

Compute matrix elements

Once the single electron orbital wavefunctions have been imported into the WFunction and WFunctions_data class objects, and the emitter configurations are defined, we are now ready to compute the matrix elements corresponding to the absorption and emission, i.e. \(v_{k,\alpha}^{(S)}\) and \(v_{k,\alpha}^{(A)}\). In practice this computation can be done in three different ways listed below with the order of increasing scalability:

(1) Computation for single multipoles

To quickly check calculation of a single case in python you can directly run the V12_par method for computation using wavefunctions, V12_ED for matrix element computation for electric dipolar transition, and V12_MD for magnetic dipolar transition. Please refer to the documentation for more details. As an example, the V12_par method can be invoked by using

v = V12_par(
    initial_wavefunction_grid: np.ndarray, # Wavefunction grid data for initial state
    k_initial_wavefunction_grid: np.ndarray, # k-operated wavefunction grid data for initial state
    final_wavefunction_grid: np.ndarray, # Wavefunction grid data for final state
    k_final_wavefunction_grid: np.ndarray, # k-operated wavefunction grid data for final state
    rgrids: np.ndarray, #Positiongrid including wavefunction grid and relative positions between multipole center and the emitter center
    k: float, # wavenumber of the photon in the medium
    L: int, # Orbital angular momentum of the photon
    M: int, # Projected total angular momentum of the photon
    P: int, # Parity of the photon
    nIndex: float = 1., # Refractive index of the medium
    Cgauge: int = 0, # Arbitrary gauge parameter. C=0 Coulomb.
    Hint: str = "Adotp", # Interaction hamiltonian type
    rshift: np.ndarray = np.zeros(3), # Additional shift of position as a control parameter
    radtype: str = "J", # Radial function of the photon mode.
    dV: np.ndarray = None, # Volume element of WF grid
    etaeff: float = 1., # Effective mass scaling
    parallel: bool = True, # execution platform- parallel or serial
    ):

In a similar way, V12_ED and V12_MD can be used when the emitter transition is known to be of ideal electric dipole or magnetic dipole types. For example,

v = V12_ED(
    pEDs: np.ndarray, # Transition dipole moment vector for the emitter transition
    rvects: np.ndarray, # Relative position vector between the emitter and the multipole center
    k: float, # wavenumber of the photon in the medium
    L: int, # Orbital angular momentum of the photon
    M: int, # Projected total angular momentum of the photon
    P: int, # Parity of the photon mode
    nIndex: float = 1., # Refractive index of the medium
    radtype: str = "J", # Radial function of the photon mode.
    Cgauge: int = 0, # Arbitrary gauge parameter. C=0 Coulomb.
    etaeff: float = 1., # Effective mass scaling
    includespin: bool = True, # Include spin in the calculation
    )

and

v = V12_MD(
    pMDs: np.ndarray, # Transition magnetic dipole moment vector for the emitter transition
    rvects: np.ndarray, # Relative position vector between the emitter and the multipole center
    k: float, # wavenumber of the photon in the medium
    L: int, # Orbital angular momentum of the photon
    M: int, # Projected total angular momentum of the photon
    P: int, # Parity of the photon mode
    nIndex: float = 1., # Refractive index of the medium
    radtype: str = "J", # Radial function of the photon mode.
    Cgauge: int = 0, # Arbitrary gauge parameter. C=0 Coulomb.
    etaeff: float = 1., # Effective mass scaling
    includespin: bool = True, # Include spin in the calculation
    )

v[i,j] represents the matrix elements for different spin flip configurations: - i,j = 0,0: spin conserving transition at both initial and final states - i,j = 0,1: spin conserving transition at initial state, spin flip at final state - i,j = 1,0: spin flip transition at initial state, spin conserving transition at final state - i,j = 1,1: spin flip transition at both initial and final states

(2) Computation for all multipoles and all positions

While the above way computes the matrix elements for a specific multipolar mode of photon, in computing the NRET rates one needs to sum over all multipolar modes, i.e. \(L = 1, 2, ... L_{max}\), \(M = -L, -L+1,..., L\) and \(P = -1 or 1\). This can be done using the computeV method. An example usage is follows:

v = computeV(
    absorber_id: str, # unique emitter-id corresponding to the absorber
    source_id: str, # unique emitter-id corresponding to the source
    absorber_transitions: list of str, # List of transitions at absorber
    ws: list of float, # List of energies of photons (in rad/s)
    Lmax: int = 2, # Max order of the multipoles
    Em_data: Emitters_data, # Emitters data object
    Pos_data: Positions_data, # Positions data object
    WF_data: list of WFunctions_data, # List of WFunctions_data objects that will be searched for the states
    parallel: bool = False,
    )

In this case the return (v) is an object of the class Vdata and specific matrix elements can be retrieved as

v_LMP = v.V[absorber_id][source_id][absorber_transition][L][M][P]

The shape of v_LMP is (2,2,) + (Nw,) + Pos_data.Rg["absorber_id"]["source_id"]. The first two axis encodes the information on whether there is a spin conserving or spin flip process at the initial and final states as also discussed under item (1) above.

(3) Compute for all multipoles and all positions in a parallel setting

To run calculations in parallel, first save the Emitters_data and Positions_data class objects containing the configurations of the emitter transitions and the spatial positions into json files by using the inbuild Encode methods. For example

with open("Em_data.json", "w") as f:
    json.dump(Em_data.Encode(), f)

with open(f"Pos_data.json", "w") as f:
    json.dump(Pos_data.Encode(), f)

Then, the Vin class can be used to create an input file corresponding to the matrix element calculation setting, for example:

vin = Vin(
    emdatafile, # path to Em_data.json
    posdatafile, # path to Pos_data.json
    wfdatafiles, # list of paths to WFdata.json files
    absorber_id, # absorber_id,
    emitter_id, # source_id,
    absorber_transitions, # list of transitions
    nIndex, # refractive index of medium
    Lmax, # Max multipolar order
    ws, # Photon energies (rad/s),
    savepath = "./",
    savefile = "V.json",
    radtype_abs = "J",
    )

vin.write_file(f"v.in")

Finally, the job can be run using the run.py script. For example,

$ srun python <path to pyRET>/run.py -in v.in >v.out

Alternatively,

$ srun python pyRET_compute_matrix_elements -in v.in >v.out

The output will be a V.json file that can be decoded to the matrix element objects by using

with open("V.json", "r") as f:
    V = Vdata().Decode(json.load(f))

References

[1] S. Chattaraj and G. Galli, Energy transfer between localized emitters in photonic cavities from first principles, Phys. Rev. Research 7, 033229 (2025).

[2] S. Chattaraj and G. Galli, Energy transfer between localized emitters in photonic cavities from first principles, Phys. Rev. Research 7, 033229 (2025).