Nonlinear Magneto-Optical Rotation with Frequency-Modulated Light

When linearly polarized light interacts with an atomic transition in the presence of a magnetic field, the polarization angle of the light can be rotated. When the rotation angle depends on the light intensity, the effect is called nonlinear magneto-optical rotation (NMOR), and is discussed in the tutorial Zeeman Structure: Nonlinear Magneto-Optical Rotation.
In a standard NMOR experiment, a zero-field magnetic resonance is observed: the optical-rotation signal is zero at zero magnetic field, increases as the magnetic field (and Larmor frequency) is increased — up to some resonance width generally determined by some relaxation rate in the system — and then falls to zero at high fields. On the other hand, if the light field is modulated, additional magnetic resonances can be observed at harmonics (and half-harmonics) of the modulation frequency. This can be very useful — for example, it allows the high sensitivity of a narrow-resonance zero-field measurement to be translated to measurements at higher fields.
In this tutorial, we model NMOR with frequency-modulated light, and show how a Floquet technique can be used to directly find the Fourier components of the density matrix in a periodic solution to the density-matrix evolution equations.
Load the package.
In[1]:=
Click for copyable input
The Floquet technique that we use below works on real variables. In order for it to work with the density matrix elements, we specify that they should be written in terms of their real and imaginary parts, using the ComplexExpandVariables option for DensityMatrix.
Specify options for DensityMatrix.
In[2]:=
Click for copyable input
Out[2]=
Specify options for plotting.
In[3]:=
Click for copyable input
We will work with a J=1→J=0 atomic transition.
Define the atomic system.
In[4]:=
Click for copyable input
Out[4]=
We want to define an optical field with a modulated frequency, i.e., where the frequency itself is a function of time. We need to be careful when doing this — the time dependence of a field with constant frequency is written as something like e-i t, and it is tempting to simply replace in this expression by our chosen time-dependent frequency (t). However, the factor t in the expression for a static field represents the accumulated phase since time zero, i.e., the integral of the frequency from time 0 to t. (Conversely, the instantaneous frequency is the derivative of the time-dependent phase.) Thus we should replace t by .
Define a linearly (x)-polarized optical field with frequency and electric-field amplitude written in terms of the Rabi frequency (up to a numerical factor) R.
In[5]:=
Click for copyable input
Out[5]=
Define a modulated frequency with central frequency , modulation frequency m, and modulation depth m, and integrate to find the total phase accumulated starting at t=0.
In[6]:=
Click for copyable input
Out[6]=
Replace t with the integrated phase in the definition of the optical field.
In[7]:=
Click for copyable input
Out[7]=
In addition to the optical field, we assume that the atoms are subject to a z-directed magnetic field with Larmor frequency L.
Define the Hamiltonian.
In[8]:=
Click for copyable input
Out[8]//MatrixForm=
Draw the level diagram, showing optical couplings and Zeeman shifts.
In[9]:=
Click for copyable input
Out[9]=
To perform the rotating-wave approximation, we want to transform into a frame that is rotating so as to acquire the same accumulated phase as found above. Since RotatingWaveApproximation takes a frequency as its third argument, we call it with an effective frequency obtained by dividing the total accumulated phase by t. This quantity will then be multiplied by t when the transform matrix is constructed.
Call RotatingWaveApproximation, which automatically generates an appropriate unitary transform using the supplied effective frequency.
In[10]:=
Click for copyable input
Out[10]//MatrixForm=
Note that in performing the transformation of the Hamiltonian to the rotating frame, the time derivative of the transform matrix is taken, which introduces a term equal to the derivative of totalPhase into the diagonal term of the excited state energy. This term is precisely the time-dependent frequency that we originally chose. In fact, the effective Hamiltonian that we have obtained is identical to what we would have gotten if we had naively performed the RWA on a Hamiltonian for a static field, and then substituted our chosen expression for instantaneous frequency in place of .
In order to generate the density-matrix evolution equations, we also need to describe the relaxation processes in the system. The excited-state spontaneous decay rate is ; we also assume a transit relaxation rate of t.
Find the matrix describing spontaneous decay and transit relaxation.
In[11]:=
Click for copyable input
Out[11]//MatrixForm=
Find the matrix describing repopulation of the ground state due to atom transit and spontaneous decay from the excited state.
In[12]:=
Click for copyable input
Out[12]//MatrixForm=
Generate the density-matrix evolution equations.
In[13]:=
Click for copyable input
Out[13]//Short=
The DM-evolution equations have explicit time dependence due to the frequency-modulation term. They can be solved numerically using NDSolve if we supply initial conditions and numerical values for all of the parameters.
Generate initial conditions corresponding to an unpolarized ground state.
In[14]:=
Click for copyable input
Out[14]=
Choose numerical values for all parameters.
In[15]:=
Click for copyable input
Solve the system of equations with NDSolve and plot the imaginary part of one density-matrix element, corresponding to an optical coherence, as a function of time.
In[16]:=
Click for copyable input
In[17]:=
Click for copyable input
Out[17]=
The solution displays an initial transient response, and then settles down to a periodic condition, with period equal to the modulation period 2/m. If we are primarily interested in this periodic state, it can be much more numerically efficient to solve for it directly, by finding the coefficients of a Fourier expansion of the density matrix, using a Floquet technique.
For the case considered here, the density-matrix-evolution equations can be written in matrix form as
where A and B are time-independent matrices, c is a time-independent vector, and here we consider (t) to be a real column vector consisting of the real and imaginary parts of the density-matrix elements.
Because the desired solution for the density matrix is periodic, we can perform a Floquet analysis of the system by expanding (t) in harmonics of the modulation frequency:
where the vectors an are the Fourier coefficients of the expansion, assumed to be time-independent for a periodic solution. (Because (t) is real, we have a-n=(an)*.) After substituting this expansion in the evolution equation, the result can be written as
Due to the orthogonality of the Fourier basis functions, this equation implies the recursion relation
where we have defined , where I is the identity matrix.
The recursion relation can be solved for the coefficients an in various ways. One method, implemented here, finds the solution in terms of a continued fraction involving the matrices B and . Choosing a cutoff harmonic, above which all Fourier coefficients are assumed to be zero, allows a numerical solution to be obtained.
EquationsToMatrixRecursion[eqs,vars,,t,n]convert a system of differential equations with explicit sinusoidal time dependence at harmonics of to matrices specifying a recursion relation for the Fourier components with harmonics labeled by n
MatrixRSolve[mats,{n,ncut,nreturn},{vars,t}]solve the recursion relation returned by EquationsToMatrixRecursion, assuming Fourier harmonics of order ncut or higher are zero, returning values for harmonics up to nreturn
FourierExpandVariables[vars,,t,nmax]write a list of variables in terms of Fourier coefficients up to harmonic nmax
FourierExpandAndIntegrate[expr,vars,,t]find the time-average of an expression over an oscillation period, and expresses the result in terms of Fourier coefficients

Solving equations using Floquet analysis and the matrix-continued-fraction method.

Convert the evolution equations to a matrix recursion expression for the Fourier coefficients.
In[18]:=
Click for copyable input
Display the form of the generated matrices.
In[19]:=
Click for copyable input
Out[19]=
Solve the recursion relation, using a cutoff harmonic of 10, and return values for the Fourier coefficients up to the ninth harmonic.
In[20]:=
Click for copyable input
Out[20]//Short=
The same DM element plotted above, written in terms of Fourier coefficients up to the ninth harmonic.
In[21]:=
Click for copyable input
Out[21]=
Compare the solutions obtained from NDSolve and MatrixRSolve. After the transients decay, the agreement is quite good.
In[22]:=
Click for copyable input
Out[22]=
Here we compare the solutions using a variable number of Fourier components.
In[23]:=
Click for copyable input
Out[23]=
In many experimental situations, lock-in detection is employed, which picks out a single harmonic from the periodic evolution. In this case, we may only need the lowest or the first few harmonics from the solution, even if higher harmonics are required to completely describe the time dependence of the system.
As an example, consider the optical-rotation signal measured with the modulated light beam. First we find the expression for the time-dependent optical rotation in terms of the density-matrix elements. We do this by calling Observables with the effective frequency totalPhase/t. We can set 0 in the resulting expression, and set m to zero, since it is very small compared to . The wavelength here contributes here only an overall multiplicative factor, and so can be set to one.
Call Observables to find the optical rotation signal.
In[24]:=
Click for copyable input
Out[24]=
We will assume that lock-in detection is used to find the in-phase and quadrature components of the signal at the first harmonic. These are found by multiplying the rotation signal by Cos[m t] and Sin[m t], respectively, and then averaging over a modulation period. The function FourierExpandAndIntegrate expresses these averages in terms of the Fourier coefficients.
Find the in-phase and quadrature components of the optical rotation signal in terms of the first-harmonic Fourier coefficients.
In[25]:=
Click for copyable input
Out[25]=
Only the first-harmonic Fourier coefficients are required in order to calculate this signal.
Here we plot, for a fixed modulation frequency, the in-phase and quadrature components of the first-harmonic signal as a function of Larmor frequency.
Choose parameters, calculate a table of the rotation signal as L is varied, and plot the results.
In[26]:=
Click for copyable input
Out[26]=
Out[26]=
The zero-field resonance analogous to the one occurring in NMOR with static light is seen at the center of the in-phase (top) plot. Additional resonances are seen at higher fields in both the in-phase and quadrature plots, occurring when twice the Larmor frequency is equal to a harmonic of the modulation frequency.
We can visualize the mechanism leading to these additional resonances by plotting the angular-momentum probability surfaces for the ground-state density matrix.
Form the ground-state density matrix by selecting levels with label 1.
In[27]:=
Click for copyable input
Out[27]//MatrixForm=
First we consider the zero-field resonance. Here we can simplify to the case of no modulation. In this case, the Hamiltonian is time-independent, and the density matrix eventually reaches a steady state. We can write the equations for the steady state condition by setting the time derivatives to zero.
Set the first derivatives of all of the density-matrix elements to zero in the evolution equations.
In[28]:=
Click for copyable input
Out[28]//Short=
These equations can be directly solved for the steady-state density matrix.
Choose parameters, solve the steady-state equations, and substitute into the density matrix.
In[29]:=
Click for copyable input
When we plot the angular-momentum probability surface, we subtract off most of the unpolarized state from the density matrix, so that the polarized part of the density matrix can be seen clearly.
Plot the angular-momentum probability surface.
In[30]:=
Click for copyable input
Out[30]=
We see that the steady state has polarization (alignment) at an angle to the x-axis. This is the result of the competing actions of optical pumping, which produces alignment along the x-axis, Larmor precession, which rotates the alignment about the z axis, and ground-state relaxation, which causes the alignment to decay after it has precessed through a small angle. This steady-state rotated alignment causes the optical rotation of the light, producing the measured signal.
As the magnetic field is increased, the rotation angle of the alignment and the light initially increases. However, for higher values of the magnetic field, the alignment can precess through 180° or more before decaying. This causes the polarization to be washed out in the steady state, resulting in a decrease in the optical rotation signal.
Calculate the steady-state density matrix for higher magnetic field, and plot the AMPS.
In[31]:=
Click for copyable input
Out[31]=
This initial linear increase and then decrease of the signal leads to the characteristic dispersive shape of the zero-field resonance (the signal flips sign if the field is reversed, because the polarization precesses in the opposite direction). Note that the steady-state rotation signal with no modulation actually has the opposite sign as the first-harmonic signal plotted above, due to specific properties of the polarization evolution.
Plot the zero-field NMOR resonance for unmodulated light.
In[32]:=
Click for copyable input
Out[32]=
For unmodulated light, the washing-out effect for high magnetic fields means that the signal is always zero at high field. Now consider the case when the light is modulated. Then the rate of optical pumping is modulated — for simplicity we can think of the atomic polarization being created in short regular pulses. If the Larmor frequency happens to equal the modulation frequency (or half of the modulation frequency, because of the two-fold symmetry of the atomic polarization), the Larmor precession will be synchronized with the modulated optical pumping. The polarization will not wash out, because the optical pumping is timed to reinforce, rather than cancel out, the existing polarization. This results in a periodic state of rotating alignment.
We can calculate and plot this periodic state using the Floquet technique employed above.
Write the ground-state density matrix as a sum of Fourier coefficients.
In[33]:=
Click for copyable input
Out[33]//Short=
Choose parameters and calculate the periodic state.
In[34]:=
Click for copyable input
Plot a series of AMPS showing the time evolution of the periodic state.
In[35]:=
Click for copyable input
Out[35]=
When the Larmor frequency is exactly in sync with the modulation frequency, the rotating alignment produces no in-phase optical rotation signal, because the atomic polarization is aligned along the optical polarization axis in phase the with the light modulation. However, there is a quadrature signal, seen in the above plots as a function of Larmor frequency. If the Larmor frequency is increased or decreased, the Larmor precession will lead or follow the modulation frequency, and an in-phase signal can be seen.
If the Larmor frequency is tuned far enough off of resonance, the polarization will again wash out and become static:
Plot the evolution when tuned away from the high-field resonance.
In[36]:=
Click for copyable input
Out[36]=
Thus the in-phase high-field resonances are dispersive, as for the low-field resonance.