Propagation of optically tunable coherent radiation in a gas of polar molecules

Coherent, optically dressed media composed of two-level molecular systems without inversion symmetry are considered as all-optically tunable sources of coherent radiation in the microwave domain. A theoretical model and a numerical toolbox are developed to confirm the main finding: the generation of low-frequency radiation, and the buildup and propagation dynamics of such low-frequency signals in a medium of polar molecules in a gas phase. The physical mechanism of the signal generation relies on the permanent dipole moment characterizing systems without inversion symmetry. The molecules are polarized with a DC electric field yielding a permanent electric dipole moment in the laboratory frame; the direction and magnitude of the moment depend on the molecular state. As the system is resonantly driven, the dipole moment oscillates at the Rabi frequency and, hence, generates microwave radiation. We demonstrate the tuning capability of the output signal frequency with the drive amplitude and detuning. We find that even though decoherence mechanisms such as spontaneous emission may damp the output field, a scenario based on pulsed illumination yields a coherent, pulsed output of tunable temporal width. Finally, we discuss experimental scenarios exploiting rotational levels of gaseous ensembles of heteronuclear diatomic molecules.


S1 Relation between broken inversion symmetry and dipole moment
A diagonal element of the dipole moment operator with i = e, g is given by where r = {r α } α=1,...,N represents the set of positions r α of all N charges q α contributing to the dipole moment of the system. For simplicity, in this proof we assume a single charge q at position r 1 . The proof for multiple charges is a straightforward generalization.
Figure S1. Spherical coordinate system used in calculations according to (S2). The unusual choice of directed angles θ and φ (notice the arrows in the picture) was intentional to set symmetric limits of integration in the equations.
To evaluate this quantity, we express r 1 in spherical coordinates defined as follows (Fig. S1): x = r cos φ cos θ , y = r sin φ cos θ , z = r sin θ . By substituting angles φ → −φ , θ → −θ in the second term and making use of the relation cos(−θ ) = cos(θ ), we obtain where −r 1 = (r 1 , −φ 1 , −θ 1 ). If the wave function r 1 |i has inversion symmetry, | r 1 |i | = | −r 1 |i | holds, and the bracket at the end disappears. This implies that the diagonal elements d ee and d gg are equal to zero. This means that broken inversion symmetry is required to create a permanent dipole moment. For the off-diagonal elements, we simply have regardless of the wave function's symmetry. Due to the odd character of the position operator, if both wave functions share the same symmetry, the integral is equal to zero, and the transition is electric-dipole-forbidden.

S2 Bloch-Maxwell equations
Here, we provide a detailed derivation of the Bloch-Maxwell equations. First, let us insert the form of the field given by Eq. (1) into the master equation (6). The set of equations for the density matrix elements has the following form: The dynamics of the other elements can be found based on the density matrix properties: ρ ee + ρ gg = 1 and r ge = r eg . Next, we make use of the ansatz (9) and expand the κ-dependent term into a series of Bessel functions e −iκ sin x = ∑ +∞ n=−∞ J n (κ)e −inx : If Ω R ω, the oscillatory terms on the right-hand side of the above equations make negligible contributions. We neglect them under the rotating wave approximation, in which out of the infinite sums the only surviving terms correspond to n = 0, 2 if they are multiplied by the drive envelope E , or correspond to n = 1 if they are multiplied by the signal field. We arrive at the equations (10) from the article. To derive formula (12), we insert Eq. (1) into the right-hand side of Eq. (4). Performing all the derivatives, we obtain Clearly, the slowly varying terms of the expression correspond to the signal, while the terms related to the drive envelope are multiplied by rapidly oscillating factors. Let us now investigate the right-hand-side of the wave equation with the second time derivative of polarization. The fastest method is to apply the derivatives to the form of polarization given by Eq. (11). We introduce the symbol The second time derivative of the polarization is The second term in the bracket reads where we have explicitly inserted the first and second derivatives of R in accordance with Eq. (S9). We insert this result back into Eq. (S10). Now, in analogy with the procedure conducted for the left-hand side (S8) of the wave equation, we can separate on the righthand-side terms proportional to different powers of the oscillating factor exp[i(kz − ωt)]. In the rotating wave approximation, we neglect powers other than 0 and 1, as only these two are present on the left-hand side given by Eq. (S8). Now, we make the important assumption that any cross-talk between the slowly varying terms and those oscillating at the frequency ω can be neglected: The polarization terms oscillating at ω are coupled to the drive, while the slowly varying terms act as the source for the signal. This assumption is valid if ∂ E signal ∂t ωE signal and similarly for the slowly varying part of the polarization. As a result, we can separate two wave equations, describing respectively the dynamics of the drive and of the signal. However, we assume the drive to be strong enough not to be affected by the coupling to the medium. The only relevant field equation is therefore the one for the drive, given by (12). On its right-hand side appear the terms corresponding to n = 1 in Eq. (S11).

S3 Numerical approach
The solution of the coupled set of Bloch-Maxwell equations (10, 12) can be found numerically.
The wave equation has the form where f is the function we wish to find, s is a source term, and c is the speed of the envelope in vacuum. Introducing discretization, we end up with a time-space grid of evenly distributed points (z j ,t i ) with respective spatial and temporal steps ∆z and ∆t. Hence, for the space variable, we have z j = z 0 + j∆z, where j ∈ [0, 1, ..., N L ], and N L is the number of the point at the end of the sample of length L, so N L ∆z = L. Similarly, t i = t 0 + i∆t for i ∈ [0, 1, ...] but with no upper bound. For convenience, we denote values of the functions at grid points f (z j ,t i ) ≡ f i, j and s(z j ,t i ) ≡ s i, j . By the solution at time t i+1 , we understand the set of values { f i+1, j } for all j, and hence we have to find an expression that depends only on the previously calculated values. This can be done by expressing the second-order derivatives by the three-point (midpoint) formula 1 . Because we solve a second-order partial differential equation, two initial conditions for each space point are required.
We introduce values f j,0 , s j,0 and velocities ∂ f j,0 /∂t ≡ g j for the initial time t 0 . To reach an accuracy on the order of (∆t) 2 , we express f 1, j as 2 and we find

3/6
Obviously, the foregoing expressions are not valid at the ends of the sample. There, we apply transparent boundary conditions 3 for the radiation exiting the sample (for j = 0, N L ): where we assumed no sources at the ends of the sample. In addition, we have introduced the parameter η = c∆t/∆z, the so-called Courant number first described in Ref. 4, whose value is crucial for numerical stability. In our case, η = 1 corresponds to the so-called "magic step" 5 and leads to a solution that is not affected by the numerical dispersion problem. The source term s(z,t) corresponds to the solution of the Bloch Eqs. (10). The latter is a set of first-order differential equations of one variable for each point in space. At each time step, we calculate the solution using the Python build-in method odeint from the scipy.integrate library. It is a RK4 method, snd so the order of accuracy is (∆t) 4 . This high level of accuracy is required because the source term is eventually differentiated twice, reducing the accuracy to the order of (∆t) 2 . The differentiation is performed using the three-point (midpoint) method.
Our solver is written in Python 2.7, where two sets of equations (S14, S15) were implemented. A comparison between Eqs. (12) and (S12) reveals that f ≡ E signal , and s is the right-hand side of Eq. (12). The solver is available on a public repository 6 . We draw the readers attention to another existing solver 7 . Figure S2. Energies of the rotational states in the ground electronic and vibrational state of the LiH molecule, as a function of electric field. The dashed line represents E DC = 150 kV/cm. Levels represented by the blue and green lines were selected for the calculations in the main text.

S4 Lithium hydride example
To perform calculations based on a real system we have chosen the LiH molecule in ground state X 1 Σ + . We investigate its rotational states |NM . Since in the ground state the electronic angular momentum is zero, N and M denote the rotational angular momentum of the molecule and its projection. The permanent electric dipole moment in the molecule's frame is d 0 = 5.88 D 8 and the rotational constant B e = 7.513 cm −19 . In general, the orientation of a molecule in the gaseous ensemble is random and as the result, in the lab frame, the average dipole moment cancels out. To distinguish one of the directions (e.g. z-axis in the lab frame) and orient molecules we may apply an additional, constant electric field E DC . This results in the DC where the spherical harmonic Y 10 has been used. The 3 − j Wigner symbols may now be used to evaluate the transition elements. As a result, the following matrix representation of the Hamiltonian is obtained where A diagonalization of the presented Hamiltonian allows us to find new eigenstates and eigenenergies of the system for different values of the DC electric field. The new eigenstates are superpositions of the original states |NM , and for the studied range of fields each eigenstate has one dominant contribution. The eigenenergies are presented in Fig. S2 for the amplitude E DC in the range between 0 and 400 kV/cm. The label corresponds to the state |NM whose contribution to the eigenstate dominates. As can be seen, the energy gap between the levels of interest grows with the electric field in the considered range.
The new eigenstates are characterized with a nonzero permanent electric dipole moment in the lab frame, oriented along the z-axis (the direction of the E DC ) as presented in Fig. S3. The black dashed line indicates the field amplitude used in the main text. Additionally, transition dipole moments between pairs of states can be induced. In Fig. S4 we present all components of the transition dipole moments between the ground state to a set of possible excited states. As expected, for M = M = 0 the transition dipole moments are parallel to the z-axis while transitions between levels with M = 0 correspond to dipoles in Figure S4. Projections on the laboratory frame axes of the transition electric dipole moments for transitions between the ground and selected excited states. The dashed line represents the amplitude E DC = 150 kV/cm. The orange curve representes the results for the pair of states selected for the calculations in the main text.
oriented in the xy plane. These dipole moments decrease with the growing field due to the increasing contribution of states different than the dominant one.