Floquet description of Optically Pumped Magnetometers

We present theoretical description of Voigt and Faraday effect based optically pumped magnetometers using the Floquet expansion. Our analysis describes the spin-operator dynamics of the first, $\hat{F}(t)$, and second, $\hat{F}^2(t)$, order moments and takes into account of different pumping profiles and decoherence effects. We find that the theoretical results are in good agreement with the experimental demonstrations over a wide range of fields and pumping conditions. Finally, the theoretical analysis presented here is generalized and can be extended to different magnetometry schemes with arbitrary pumping profiles and multiple radio-frequency fields.

Most of OPMs are based on a Faraday dispersive measurement, in which oriented states (see Fig. 1 (a)) are prepared and probed by a detuned laser beam measuring the Faraday rotation induced by the spin polarized sample. This kind of configuration can run in scalar or vector mode [17,18], typically in orthogonal geometry. A different approach has been shown in ref. [19], in which an aligned state is prepared instead of an oriented one (see for example Fig. 1 (b)) and read through paramagnetic resonance i.e. non-dispersive measurement. This kind of state allows a vector magnetometer operation using radio-frequency fields [20], or, as it was proposed more recently [21], adopting an all-optical approach which performs a dual axis magnetometer based on Hanle effect. On the other hand, in ref. [22] we have shown that indeed it is possible to employ dispersive measurements based on Voigt rotation when working with aligned states driven by radio frequency fields, also showing vector magnetometry operation, (see Fig. 2). However, the typical description for oriented spins probed by the Faraday rotation, which depends on the first moment of the spin operator is not suitable for the aligned states probing the Voigt * hans@if.usp.br † tadas.pyragius@nottingham.ac.uk rotation as the latter is proportional to second moment spin operator. Another important and common feature of the oriented and aligned spin state magnetometers is the use of synchronous pumping in which the amplitude (or frequency) modulation avoids optical decoherence due to the pump. This implies that the pumping rate and the decoherence rate (e.g. square waves intensity profile) are in general time dependent. This type of time modulation leads to interaction with many harmonics, which is typically avoided, by assuming that the pumping is weak, or approximating the interaction by neglecting higher harmonic terms. In this work, we show an approach that is capable of encompassing these previously neglected time dependent terms. To do so, we employ a Floquet expansion to solve the spin dynamics which enables us to build a more realistic picture of the spin evolution [23,24]. The spectral decomposition of the Floquet expansion allow us to address the solutions for multiple harmonics generated in the dynamics independently, which can be directly compared with the experimental observations.
In this paper we present the dynamics of the first and second moment elements when the spins are driven by a radio frequency field, which describes the dynamics of oriented and aligned states. Furthermore, we find the solution for the spin dynamics in the realistic situation arXiv:2011.03785v1 [quant-ph] 7 Nov 2020 FIG. 2. The laser interacts with the atomic ensemble allowing a dispersive measurement of the aligned state dynamics. The dynamical evolution of aligned states dressed by a radiofrequency field enables detection of all three vector components of the magnetic fields. The magnetic fields Bx, By and Bz are measured from the change in ellipticity on the probe beam and can be represented as a frame rotation from (x , y , z ) to (x, y, z) [22].
where the optical pumping presents an arbitrary time dependence. We show that the general dynamics of the second moment in the Liouville space can be reduced to Bloch equations (see eq. (2)) and can be solved by employing the Floquet expansion, which to our knowledge has not been reported before. This solution predicts sensitivity to all three vector components of the magnetic field as reported in the experimental work in [22]. This approach is also compatible with scenarios in which multiple radio-frequency fields are used.
The paper is organized as follows. In Section II we introduce the differences between the Faraday and Voigt rotation in the dispersive regime with respect to the statistical moments of the spin operators. In Section III, we study the spin dynamics including the interaction with static and radio-frequency magnetic fields, and consider an arbitrary time dependence for the optical pumping. Section IV describes the Floquet expansion to solve the dynamics of the first moment. Section V describes the dynamics of the second moment of the spins and the transformation to the Liouville space. Section VI shows the results on the Floquet expansion to solve the dynamics and we discuss the main features of the model. Section VIII presents our conclusions.

II. DISPERSIVE OPTICAL MEASUREMENTS
The most common configuration for linear optical devices uses dispersive optical measurements of spinpolarized atoms to detect the presence of an external magnetic field. The modulation of the birefringence of the medium caused by the Larmor precession of such spin-polarized atoms can be observed polarimetrically, this is known as Faraday rotation. More specifically, in terms of polarization moments, the spin oriented atoms represented by the probability surface in Fig. 1 (a), for which the z-component of the total angular momentum F = (F x ,F y ,F z ) on average satisfies F z (t) = 0, induces a polarization rotation of an incident linearly polarized light that propagates along the z-direction [25,26]. Assuming that such interaction is with an atomically thin sample with no back action effects, the rotation can described by the Stokes' parameter. For linearly polarized such light-matter interaction is given by where the prime indicates the output field after interacting with the atomic medium, represent the photon flux of elliptical and at 45°polarized light expressed in terms of creation and annihilation operatorsâ ± andâ † ± for circular polarization components; G (k) F is the rank-k coupling strength and n F are the atoms with the same F -manifold state [27].
In terms of statistical definitions, the average value of the spin operator F i (t) where i = x, y and z, correspond to the first moment of a statistical distribution of the spin operatorF i (t). Thus, we can claim that OPMs based on the dispersive Faraday rotation work only for quantum states with non-zero first moments. Since F i (t) correspond to the classical description of polarized samples, the dynamics of OPMs based on rf excitation or optical excitation, are classically described by the Bloch equations [28] where γ is proportional to the Larmor frequency and the polarization vector is P = (P x , P y , P z ) with P i = F i . The decoherence T 1 accounts for relaxation of the spins along the longitudinal direction from magnetic field gradients, collisions with the walls and atoms that are pumped with spin polarization P 0 . The decoherence of the transverse polarization is described by a term proportional to 1/T 2 , and represents the atom-atom collisions. When the spins are driven by a sinusoidal magnetic field B the spin dynamics follows a resonance response, which is given by the Bloch solution [28,29]. We can also consider the precession of aligned states as the one represented in Fig. 1 (b), around a static magnetic field. To optically probe this kind of dynamics dispersively, we have proposed in ref. [22] a measurement based on Voigt rotation. This effect measures the changes in the linear birefringence of the probing light (see Fig. 2). In the limit of a far detuned probe light traversing an atomically thin sample and assuming no back-action effects, the changes in the linear birefringence can be described by the following Stokes' parameters [30] Ŝ which is proportional to transverse second moments of the spin operators i.e the average of the second order products F i (t)F j (t) with i, j = x, y and z. Unfortunately, the dynamics for the second moments cannot be described by the Bloch equations, eq. (2). The dynamics for the first moment P i = F i , which is the average of a linear operator, in general is different from the dynamics for the second moments F i (t)F j (t) , which is the average of bi-linear operators. In particular, there is no Faraday rotation for an aligned state out of the orientedto-alignment conversion (OAC) regime, which yields a trivial solution for the Bloch equation. Hence, it would be desirable to have a dynamical equation like eq. (2), but for the second moment operators based on Heisenberg equations of motion.

III. SPIN DYNAMICS IN A RADIO-FREQUENCY DRESSED FIELD WITH OPTICAL PUMPING AND RELAXATION
A. Heisenberg-Langevin equations Consider a radio frequency dressing field in the presence of transverse and longitudinal magnetic fields where B rf is the amplitude of the rf field and B dc is the static field along the longitudinal direction, which are experimentally controlled. Additionally, we have the external fields B ext i with i = x, y and z, which originate from external sources. For small fields, where the second order Zeeman shift can be neglected, the magnetic field interaction is given byĤ = (µ B g F / )F · B, such that in the Heisenberg picture we havê with g F = g F / and Ω i = µ B g F B i with i = rf, dc, x, y and z. The coherent part of the atomic spin dynamics is given by the Heisenberg equation such that the spin dynamics due to the magnetic fields is given by where the atomic spin vector isF(t) = (F x (t),F y (t),F z (t)) and Now, to describe the full dynamics of the magnetometer, we need to include additional terms, which govern the pumping and the decay processes of the prepared spin states. One of the contributions corresponding to the state preparation process is the pumping term, namely where the pump rate Γ p (t) describes a general form in time at which the state preparation is done e.g. synchronous pumping with any harmonic profile. We will later describe the harmonic decomposition of Γ p (t). We consider the action of the pump process as stochastic flips in time on the atomic operator. Therefore, in general, we considerF in (t) as a stochastic vector operator with non-zero mean value. Thus, we propose a linearized-like version of the input operator i.e.Ô = O + δÔ, such thatF with non-zero mean value, where the stochastic part satisfies F in (t) = 0 and its correlation function is where σ in is the input second moment matrix, which acts as a diffusion term in the second moment dynamics. This input operator that we propose recovers the pumping term in the Bloch eq. (2), since F in = P 0 e z . By adopting this linearized version of the input operator, which follows a perturbative approach of the pump action into spin dynamics, splits the contributions given by the first moment and the second moments of the input state. The second term we want to include is related to the relaxation of the spins. We apply a relaxation process in terms of stochastic operatorsF (t) = (F x (t),F y (t),F z (t)), to obtain the Langevin dynamics of the atomic spin operators where F (t) = 0 and in the general case the relaxation matrix Γ rel is a diagonal matrix with components Γ i with i = x, y and z. The stochastic operator associated to this relaxation process satisfies the correlation function whereΓ correspond to the diffusion matrix. In Appendix B we formulate the diffusion matrix in terms of the relaxation matrix Γ rel satisfying operator commutation relations.
In general, each direction is subjected to different decoherence rates. Nevertheless, the relaxation processes in the transverse directions are typically different from the longitudinal direction. As a result, it is commonly considered that the transverse direction is affected equally by spin exchange collisions such that Γ x = Γ y = Γ 2 . On the other hand, processes like wall collision, decoherence induced by the pump and magnetic field gradients, may relax the longitudinal direction at a different rate, given by Γ z = Γ 1 .
Combining the terms containing coherent spin dynamics, pumping and relaxation, we obtain the total spin dynamics of the system IV. DYNAMICS OF THE FIRST MOMENT F (t)

A. Spin dynamics in the laboratory frame
The dynamics of the first moment is defined by the mean value of eq. (14), such that P(t) = F (t) , which corresponds to the classical description of magnetic spins. Therefore, its dynamics can be written as where we have defined ext − Γ rel and P in = F in . The matrix B(t) can be decomposed spectrally in which Another term that can be spectrally decomposed is the the pumping rate such that, for instance, a square-wave pumping profile can be decomposed to where d corresponds to the duty cycle of the carrier wave. The general description in eq. (19) can simulate a broad range of time dependent pumping rates with different spectral decompositions e.g. sine, sawtooth etc. With the definitions above, the dynamical equation can be written as To check the consistency of this solution, we show in Appendix A that applying the rotating frame transformation and considering no external magnetic fields, we recover the Bloch solution [28]. Nevertheless, it is worth noting that eq. (21) and its counterpart in the rotating frame in eq. (A22) stands for a more general description, which accounts for the presence of external fields and more realistic pumping schemes, which could be harmonically decomposed, as a result, it requires Floquet expansion to solve it. Given the harmonic nature of this dynamical equation, we employ a Floquet expansion of the spin operatorsF(t) in order to find a steady state solution for all the possible harmonics. Therefore, we expand the spin operator harmonically aŝ such that for the first moment we have P(t) = n P (n) (t)e inωt where P (n) (t) = F (n) (t) . This expansion allows to compute the dynamics of P(t) by finding the time evolution of the harmonic components P (n) (t). Substituting this expansion into eq. (21), we find the following dynamics for the spin harmonics where B (0) n = B (0) − inωI and Q ∈ Z corresponds to the cut-off frequency index i.e. a finite, but large number of harmonics required to satisfy convergence in numerical calculations [23,31]. It is worth noting that, in the laboratory frame, the matrices B (±1) are directly proportional to the amplitude of the rf field, responsible for coupling the harmonics of P (n) and P (n∓1) , respectively. However, in Appendix A 1 we show that in the rotating frame, the harmonics are coupled by the matrices M (±1) , which are proportional to the external transverse fields. More specifically, the in-phase component of M (±1) is proportional to magnetic fields in the x-direction, whereas its out of phase component is proportional to the external magnetic field in the y-direction. This is going to be described in more detail in the results section.
Defining the vector V with V i = 1 for all i's, the spin operator can be written as P(t) = V · P(t). Therefore, the dynamics of the harmonics in eq. (23) can be written as where B = B − iNω. The calculation of the dynamics requires a limitation of this vector by defining a cut-off harmonic n = Q such that the solution converges. In the finite case, we have where the pump matrix elements p I 3×3 , and the input vector (P in ) n = P in . This new vector has dimensions d P = 3 × (2Q + 1), and the matrices have From eq. (25), the steady state solution is The convergence of this solution is verified when P F reaches stability by comparing the calculation with B 2Q+1 and B 2Q+3 , with Q being the cut-off frequency. The eq. (25) is the general form of the classical solution in eq. (2), where an arbitrary pump intensity profile is considered. In this work, we aim to calculate a more general solution for the spin dynamics, not only showing the first moment solution, but going up to second order moments of the spin dynamics. In particular, we are interested in the solution of the spin co-variance matrix when the atoms are in the presence of an external magnetic field where the spin state is prepared using a synchronous pumping process with a square wave intensity profile. This corresponds to our Voigt effect based 3D vector magnetometer described in ref. [22].
Before we discuss the first and second moment solutions of the spin operators, it is worth briefly showing the analogy between the dynamics of the first two statistical moments. We have already shown that in the laboratory frame the first moment follows the dynamics given by whilst the second moment will follow, in the Liouville space, an equivalent dynamics given by in which X(t) corresponds to the vector representation of the second moment matrix σ(t) = F (t)F(t) T in the Liouville space. This analogy allows us to easily show that the Floquet expansion employed for the first moment solution can be extended to the second moment solution and that the iterative formulas in both cases will be equivalent.

V. DYNAMICS OF THE SECOND MOMENT F (t)F(t) T
A. Dynamics of the second moments in the laboratory frame Now let us draw our attention to determine the dynamics of the second moment. To do so, we define the second moment matrix From the spin dynamics we can determine the dynamics of the second moment and substituting eq. (14) and its transpose into eq. (31), we obtain The first two terms on the right hand side of eq. (32) contain the coherent part generated by the interaction between the atomic spin and the magnetic fields. The two terms proportional to the pump rate Γ p (t) describe the dynamics of the second moment due to the pumping process. The last two terms correspond to the contribution of the stochastic noise from the pumping process and unpolarized atoms.
In particular, we need to determine the cross correlation of the atomic spinF(t) with the input and unpolarized stochastic operatorsF in (t) andF (t). In Appendix B 1, we show in eq. (B26), that the input stochastic term satisfies the following expression and the unpolarized stochastic term is given by eq. (B29) such that, substituting the equations above into eq. (32), we finally obtain In particular, for polarized samples prepared in an aligned state, F in = 0, we have B. Second moment matrix dynamics in the Liouville space The dynamical equation for the second moment is linearly equivalent to the equation for the first moment. To explicitly show that equivalence we transform the second moment matrix σ(t) into a vector X(t) in the Liouville space such that σ(t) → X(t). Meanwhile, matrices that operate in the Euclidean space from the left and from the right with respect to an operatorÔ, i.e.Ôσ A and σ AÔ , respectively, are mapped into matrices in the Liouville space such thatÔ Hence, the dynamical eq. (36) can be written in the Liouville space as where From now on we refer to X(t) as the second moment vector in the Liouville space. We have demonstrated the equivalence between the dynamics of the first and second moment that we pointed out at the end of Sec. IV B.
Since the magnetic field interaction matrix B(t) decomposes as in eq. (16), in the Liouville space this is equivalently expressed as where C (n) = L(B (n) ) + R(B (n)T ) with n = −1, 0, 1. Again, this can be solved by employing the Floquet expansion.

VI. FLOQUET EXPANSION OF THE SECOND MOMENT OF THE SPIN OPERATOR
To obtain the Floquet expansion of the dynamical eq. (39), we must check how the time dependent variables are harmonically expanded. Let us start with the harmonic expansion of the second moment matrix, which follows the expansion in eq. (22), such that in which the harmonic component σ (n) (t) can be generally expressed in terms of the average spin operator products as whereσ (n,m) (t) = F (n) (t)F (m) (t) T . Therefore, in the Liouville space, where the matrix σ(t) is represented by the vector X(t), the second moment vector is expanded as Similarly to the spin polarization P, we can now define the vectors in the matrix form for the spectral space The Floquet expansion for the dynamical eq. (39) is obtain by substituting the spectral expansions of the second moment vector X(t) and the matrices C(t) and Γ p (t). Therefore, according to eqs. (47), (19) and (42), and associating terms at the same harmonic frequency, we obtain the following recursive equation Notice that the first three terms on the right hand side are equivalent to those obtained for P (n) for the first order moments. The spectral convergence of the last term for a given n th -harmonic depends on the decay of the harmonic decomposition of the pumping rate e.g. γ (n) l ∼ Γ (n) ∼ 1/n for a square wave. As in the case of spins P, eq. (48) can be expressed in the matrix form as whereC = C − iNω takes the same form of B in eq. (26), which it components arẽ where the pump matrix term is ( Therefore the steady state is governed by From this solution we can notice that for a dominant pumping rate, in which Γ G , Λ rel , the steady state is directly proportional to the input of the second moment vector X F ≈ Γ −1 Γ in X in , losing any resonant response to the magnetic fields. In the case where collisional processes dominate the dynamics, the solution corresponds to the one for an unpolarized state X F ≈G −1 Λ rel X 0 , which also has no resonant response to the magnetic fields. Furthermore, in the case of the laboratory frame, the first harmonic components C (±1) , describe the circular components of the rf field. This does not mean that the steady state solution will contain spin oscillations at the first harmonic when the static field is on resonance. However, approaching the dynamics from the rotating frame we can determine, for some parameters, which element of the dynamics leads to some specific harmonics at resonance. In Sec. VI B we will discuss the dynamics in the rotating frame.
A. OPM response during the probe cycle From the steady state in eq. (51) one can notice that the pump rate matrix broadens the magnetic response by the Lorentzian term [C − Γ] −1 . To avoid broadening, we utilise pump-probe strategy as it is shown in Fig. 3, where the pumping is followed by a free induction decay of the spin evolution, which is where the state is probed.
FIG. 3. Free induction decay sequence of Voigt effect magnetometer. The first part of the sequence consists of state preparation (pump cycle) to prepare the aligned states. After the pump cycle a probing pulse interrogates the free induction decay dynamics. This configuration is defined as a doublestep measurement. On the other hand, probing during the state preparation is defined as a single-step measurement.
Therefore, after a pumping cycle, we have a cycle of atomic dynamics without pumping, such that integrating this expression yields the time evolution To determine the amplitude of the harmonics at the beginning of the probe cycle we integrate over a cycle of the rf frequency, such that the real and imaginary parts for a given harmonic n are where T = 2π/ω. Substituting the expansion X(t) = Q n=−Q e iωnt (X F (t)) n for the real part we have Since (X F (t )) m are slow varying envelopes of the harmonics within a time T , they are even, and therefore, the only non-zero values are Following the same procedure, the imaginary part is As a result, we need the integration of the elements X (n) (t ). From the solution in eq. (53) and considering that the initial second moment matrix is given by the steady state solution X F (0) = X s F in eq. (51), we obtain From this solution one can determine the specific harmonics that contribute to the change in ellipticity in eq. (3). Since the Floquet expansion is applied to the the atomic term F 2 x −F 2 y , the ellipticity in eq. (3) can be equally expanded as where the harmonic components h (n) X,Y correspond to the quadratures that can be measured using lock-in detection and defined in terms of eqs. (57) and (58) as in which [X I ] 1 correspond to the real and imaginary components of the nth-harmonic and the subscript 1 represents F 2 x , whereas the subscript 5 represents F 2 y . In particular, we are interested in the quadratures of the first harmonics h X , because those are the signals that can map the three components of the external magnetic field [22].
Notice that the the steady state solution in eq. (51) and the dynamical solution in eq. (53), with the use of the harmonic space can be used to construct a broad range of pumping and probing schemes (single step or double step measurement), as well as include various forms of external magnetic field interactions. This shows how the algebraic complexity of the coupling of the harmonics seen from the Liouville space, is effectively cleared in the harmonic space, allowing us to distinguish the contribution of different elements in the atomic spin dynamics.

B. Second moment dynamics in the rotating Frame
Although all the numerical results that are going to be presented in the following sections are calculated in the laboratory frame (due to its easier computational implementation), it is nevertheless worth exploring the analytical solution in the rotating frame in order to develop a physical insight into the magnetometer response.
The main feature of this vector-OPM is that the mode quadratures of the first harmonic map the transverse fields. In order to show this, let us consider the transformation in eq. (A1) applied to the second moment matrix in the laboratory frame in eq. (30) The corresponding dynamics of the second moment matrix in the rotating frame is described in Appendix C. In the Liouville space the dynamics of the second moment are given by eq. (C12) where G(t) describes the resonant and external fields, Λ rel accounts for collisional relaxation rates and pump relaxation rate Γ p (t) in the rotating frame defined in eq. (C11) Similarly to what was shown in Sec. V and VI, the G(t) decomposes as G(t) = G (0) + G (1) e iωt + G (−1) e −iωt , where G (n) = L(M (n) ) + R(M (n)T ) with n = −1, 0, 1, and in particular, G (±1) is only dependent on the transverse external fields according to the definition of M (±1) in eq. (A35). Therefore, the Floquet expansion for the dynamics of X (t) is given by whilst the pump matrix term is (Γ in ) nm = 2δ nm Γ p (n) , the relaxation matrix (Γ) nm = 2Γ (n−m) p I 9×9 and the input second moment vector (X in ) n = X in . The pump elements Γ p (n) are defined in Appendix E, eq. (C11). Moreover, the unpolarized drift matrix is (Λ rel ) nm = Λ rel for n = m = 0 otherwise is zero. Hence, the steady state solution in the rotating frame is which takes the same form as the steady state in the laboratory frame. Transforming back to the laboratory frame (see Appendix C 1), from eq. (C14), the solution is where R (n) are Liouville harmonic rotation matrices. For a general case it is not straight forward to specify the contribution of the processes into the detection of atomic response at a particular harmonic frequency. Nevertheless, we can discuss different scenarios where we can clearly determine the appearance of harmonics.
First, notice that in the rotating frame, the matrices G (±1) couple the harmonics among them, which are directly mapped onto the M (±1) matrices. According to eq. (A35), these matrices depend directly on the transverse fields, in contrast to the lab frame, in which the M (±1) matrices are only dependent on the rf amplitude. The fact of M (±1) matrices map directly the presence of the transverse fields is consistent with the experimental observation, since only transverse field can give rise to the first harmonic.
So, let us consider the case where there are no transverse fields and we have a constant pump in the system i.e. G (±1) = 0 and Γ (n) p = Γ P δ n,0 . This condition decouples the harmonics in eq. (65), diagonalizes the matrixG and therefore, the steady state solution is for the aligned input state along the rf-axis. Since the only non-zero solution is X (0) , according to eq. (68), the only harmonics present in the system are It can be shown that the transverse components of X (±1) are zero due to matrices R (±1) . Therefore, the birefringence from the atoms is described only by X (0) and X (±2) . However, when weak transverse fields are present such that G (±1) take non-zero values, the next contribution to the solution in eq. (68) are X (±1) and then which represents the presence of the weak signal at ω. This is how, for low external fields, the quadratures at ω directly map the external transverse fields.

VII. RESULTS
Here we aim to analyze the second moments of state in the hyperfine level F = 2 of alkali atoms like Rubidium which we can relate to experimental observations in our previous work [22].For this hyperfine level the cutoff frequency index is Q = 5 = 2F + 1 for which the computation converges. Nevertheless, the model can be applied to other atomic species with different hyperfine level structure.

A. Second moments for aligned states
Let us first analyze the dynamics involving synchronous pumping, in particular, for a stretched state along the x-direction, which corresponds to a mixture of equally populated states |F = 2, m F = ±2 x with its second moment vector given by X in = (4, 0, 0, 0, 1, 0, 0, 0, 1) T . Figure 4 shows the harmonics h X as a function of the static field with its Larmor frequency Ω dc . The curves describe the second moment dynamics after the state preparation in a double-step measurement as it is shown in Fig. 3. The insets show the dynamics of the observables probed in the steady state during the state preparation, which corresponds to a simultaneous pump and probe measurement (single-step measurement). Figure 4 compares the results from eq. (59) for three different situations: (a) and (b) with synchronous pumping at d = 10% duty cycle with and without the transverse external fields; (c) with continuous-wave (cw) and transverse external fields.
In the case of no external fields in Fig. 4 (a), one can notice that the real and imaginary part of the first harmonic are zero, whereas the real part of the second harmonic h (2) X has a resonant profile, reaching its maximum when the Larmor frequency Ω dc = ω rf . This is consistent with the results in sec. VI B, eq. (69) where the absence of transverse fields gives rise to rf signals only in the second harmonic. This can be described in terms of the state probability surface as shown in Fig. 5 (a). In this case, the stretched state precesses around the static field which is perfectly aligned with the quantization axis. After half a cycle, the surface returns to its initial position, in which F 2 x − F 2 y > 0, hence oscillating at twice the Larmor frequency.
In the case of a pump beam with a time dependent square intensity profile which is synchronously modulated as shown in Fig. 4 (b), with the static field being resonant with the rf driving frequency, the quadratures of the first harmonic follow a dispersive profile due to the presence of weak transverse fields, whereas the second harmonic profile shows resonant behavior, describing very closely what was observed experimentally [22]. From the rotating frame description, eq. (74) shows that for weak transverse fields the first harmonic becomes non-zero. Fig. 5 (b) shows that in this case the weak transverse fields slightly tilts the static fields axis inducing a precession of the stretched state. In this situation the second harmonic dominates, since the projection F 2 x − F 2 y > 0 oscillates at twice the Larmor frequency, but the axis of the aligned state only returns to its initial position after an entire Larmor cycle.
In contrast to the modulated pump, in the case of cw pumping, we have Γ (0) = 0 and Γ (n) = 0 for n > 0. Fig. 4 (c) shows that the first harmonic quadratures no longer display dispersive characteristics and the amplitude of the second harmonic is reduced. This means that spin evolution is not sensitive to the magnetic fields in all three directions. Now we turn to the analysis a second example. Another type of an aligned state is a pure state given by  Fig. 4 (b), except with an opposite sign. The inset shows the mode amplitude for no external magnetic field, in which the only non-zero rf signal is the second harmonic. Similar to the stretched state in Fig. 5, the aligned state |2, 0 x precesses around the static field maximizing the F 2 x − F 2 y which oscillates only at twice the Larmor frequency (see Fig. 6 (c)). Once the external transverse field is present, the axis of symmetry is tilted, and its axis precession contributes to rf-signal at the Larmor frequency. It is worth noting that the in-phase top view of the aligned state (insets) is 90 degrees rotated with respect to the stretched state in Fig 6 (c), which shows that F 2 x − F 2 y < 0 and therefore the mode amplitudes h (n) X,Y have the opposite sign.

B. Second moments for oriented states
Another effective state that maximizes the precession of the second moments is a transverse oriented state. In particular, consider |F = 2, m F = 2 x where its second moment vector is given by X in = (4, 0, 0, 0, 1, −i, 0, i, 1) T , very similar to the one for the stretched state. Fig. 7 shows the first and second harmonic as a function of the static field for an oriented state along the x-axis synchronously pumped, which achieves the same oscillation amplitude as an aligned state. The dispersive profile of X,Y is non-zero due to transverse fields, whereas the inset shows zero first harmonic response with no transverse field, as in the case of both aligned states discussed above. Fig. 7 (c) and (d) show the representation of the oriented state precessing with and without a transverse field, in the same manner as the aligned state. The inphase top view of the aligned state (inset) is not rotated with respect to the stretched state in Fig. 5, which shows that F 2 x − F 2 y > 0, and the signs of the signal are not inverted with respect to the stretched state.
The examples considered here with these three states show the atomic spin interaction with the rf field and the external field. One thing to notice is that for optimum sensitivity the aligned states are generated with a parallel pump-probe configuration, whilst the oriented stated requires an orthogonal pump-probe configuration for its optimal preparation. Hence, the aligned state are more convenient for developing miniature sensors, since they can be operated in parallel single-axis configuration.

C. Optimization
To optimize the vector-OPM response to external fields, one of the parameters to characterize is the amplitude of harmonics h (n) X,Y with respect to the duty cycle of the pump beam. Fig. 8 (a) shows the result for the normalized figure of merit (FOM)[? ] between the mode amplitude A and line-width, Γ, of the harmonics as a function of the duty cycle for the three different states analyzed above. Clearly, the three states show the same behavior of the FOM with respect the duty cycle, considering the rf amplitude as Ω rf = 0.05 ω rf . One can notice that below d = 10% duty cycle the FOM tends to zero. At d = 10%, the three signals reach an optimum point and for higher duty cycles the two curves show that the FOM decrease. It is worth noting that changing the duty cycle, without changing the pulse amplitude, reproduces the conditions the typical experimental conditions where the acousto-optical modulator (AOM) is used to generate synchronous pulses. In both cases, the changing duty cycle not only changes the duration of the pump interaction with the atomic spins, but also the effective power. Hence, for d = 0% there is no state being prepared resulting in zero Voigt rotation since the atoms are in a thermal state which has zero birefringence. In the intermediate region of the duty cycle, 0 < d < 10%, there is not enough time to prepare the optimal state for it to interact with the static fields and so the precession is small resulting in low birefringence. Around d = 10% the spin dynamics is at its optimum where the interplay between the input state and the driven rf-field result in a strong precession around the static field. Lastly, when d > 10% the pump dominates the interaction with the spin inhibiting the spin precession around the static field.
A different situation is observed with a lower rf amplitude Ω rf = 0.03 ω rf . Fig. 8 (b) shows the same behavior of the three states, but in this case the first harmonic h (1) X,Y finds an optimum FOM at d ∼ 6%, different from the second harmonic with its maximum is still at d = 10%. For higher duty cycles the FOM decays faster compared to the case with Ω rf = 0.05 ω rf . Notice that the first harmonic sensitivity to transverse fields implies that the state precesses around a tilted static field going off in the x, y plane. Therefore a slight difference in the optimum duty cycle between h (1) and h (2) relies on the fact that a weaker rf driving field requires a shorter pump time for the precession to occur going off in the x-y plane. Otherwise the pump dominates, maintaining the spin precession close to the x, y plane, which results in a high amplitude of h (2) , but a reduction of the sensitivity for h (1) .
A second parameter that determines the optimisation of the vector-OPM, is the rf amplitude. Fig. 9 (a) shows the FOM as a function of the rf amplitude of the resonance profile of the first and second harmonic. This graph describes the behavior of the three states considered in the previous examples. A high FOM for the second harmonic h (2) X is reached for very weak rf amplitudes, which corresponds to the narrowest magnetic resonance. By increasing the rf amplitude, the Voigt rotation signal at 2ω increases at the expense of broadening the magnetic resonance, which reduces the overall OPM sensitivity. Meanwhile, the maximum FOM for each quadrature h reaches its maximum at Ω rf ∼ 0.05 Ω dc . For lower rfamplitudes (below Ω rf < 0.045 ω rf ), we observe the highest relative difference between the quadratures, since the radio frequency field can not drive the tilted precession of the aligned state, whereas for higher amplitudes (above  9. (a) FOM for ω and 2ω harmonics as a function of the radio-frequency dressing field amplitude Ω rf . The FOM for h (1) X and h (1) Y is normalized by the maximum of both signals, whereas the FOM for h (2) X is normalized by its maximum on resonance. (b) Normalized mode amplitude of the second harmonic h (2) X as a function of phase of the pump beam relative to the radio-frequency driving field. Here the pump is at 10% duty cycle. The parameters used for the calculation are the same as those in Fig. 4. Here, the atoms are dressed with a uniform 5 kHz rf field and pumped with a square intensity profile with a 10% duty cycle, as in ref. [22]. The insets show the OPM response for small external fields, which behave linearly in the hx − hy − hz space. The parameters used for the calculation are the same as those in Fig. 4.
Ω rf > 0.075 ω rf ), the driving field broadens the resonance reducing the sensitivity [? ]. It can be observed that the optimization of the 3D operation is within the range near the maximum Ω rf ∼ 0.05 Ω dc in which the three components are simultaneously sensitive to magnetic fields. Within this range (light red square in Fig. 9 (a)) the three signals are non-zero and the relative FOM between the h (1) X,Y do not exceed 80%. This results in a comparable sensitivity between the external transverse and longitudinal fields.
In addition to the duty cycle and the rf amplitude optimization, Fig. 8 (b) shows the normalised amplitude of the second harmonic h (2) X a function of the pump phase with respect to the rf field. Notice the maximum amplitude of the second harmonic occurs in phase (0, 2π, · · · ), whereas a reduced and sign inverted amplitude is reached for anti-phase (π, 3π, · · · ). This curve describes the same behavior for the aligned and oriented states. The relative phase between the pump and the rf field strongly affects the state preparation process and the consequent OPM response to external fields.

D. 3D Vector mapping
In this section we show how the theoretical model predicts the three-dimensional vector field mapping operation of the magnetometer. Following the same procedure as in ref. [22], we can determine the three components of the field by setting the static field B z at B 3D z = B res + B rf /2, which maximizes the mode amplitudes h  (2) X quadratures, we are able to map the magnetometer response. In order to have a better picture of the spatial distribution of the external fields, we adopt a notation that relates the harmonics with the spatial directions x, y and z as (75) Figure 10 shows the vector magnetometer operation visualized on a 3D plot for an stretched state |2, ±2 x . Every oviform surface corresponds to the three different external longitudinal fields, B ext z . The Floquet expansion not only reproduces the oviform profile considering a large range of transverse fields as in ref. [22], but also demonstrates the vector operation for small transverse fields with the three planes in the h x − h y − h z space, which in the small field regime correspond to linear external field mapping One can observe in Fig. 10 (a) that with no transverse fields present, the spins are correctly mapped at h x = h y = 0. In this situation the aligned state precesses around the static field applied, perfectly aligned with the quantization axis of the of the probe beam, as discussed in sec. VII A. In the weak transverse field regime, the presence of orthogonal transverse fields B x and B y translates into the response of h x and h y quadratures as it is shown in insets of Figs. 10 (a) and (b). At higher transverse fields, the resonant response is broader, inducing fast transitions among the Zeeman levels, leading to a thermal state for which the second harmonic is drastically reduced in which the three surfaces tend to h z = 0, e.g. h (2) X ≈ 0 (see Fig. 10 (b)). The same kind of vector-magnetometer response describes what is observed for the aligned state state |F = 2, m F = 0 x and the oriented state |F = 2, m F = 2 x . The graphs are not shown independently to avoid repetition.

VIII. CONCLUSIONS
We have presented a theoretical model to describe the dynamics of a new kind of radio-frequency dressed three dimensional vector magnetometer based on the Voigt effect. As shown, our model describes the spin dynamics not only of the first moment, which is in agreement with the Bloch solution for Faraday based magnetometers, but also for the second moment. We demonstrated that oriented and aligned states would present vector magnetometer response by dispersive Voigt rotation measurements, however, the aligned state is compatible with parallel geometry which is more suitable for miniature sensors. In addition, we have shown that the time dependent dynamics involving synchronous pumping for state preparation can be solved employing a Floquet expansion, and the results are in agreement with the experimental observations in ref. [22]. The Floquet expansion is a powerful tool because it can solve different kinds of time dependent profiles making our approach general towards applications in understanding alternative approaches in OPM architecture. To determine the noise properties and sensitivity limits of the Voigt effect vector-OPM it would be necessary to solve the dynamics of the fourth moment of the spins. The model proposed in this work paves the way towards finding the solutions for higher moments. The time evolution of the spin operator is easily solved in a rotating frame, which oscillates at the same frequency as the rf field. In the rotating frame, the atomic state is transformed as |ψ(t) = U (t) |ψ(t) whilst the spin operators followF i (t) = U (t) −1F i (t)U (t). Considering the unitary transformationÛ (t) = e iωtFz/ ,the matrix representation of the rotating frame transformation around the z axis is written aŝ

IX. ACKNOWLEDGEMENTS
where and its inverse From eq. (A1), we can obtain the dynamical equation in the rotating frame Considering the fact that the rotation matrix R(t) can be expressed as where we have defined (A7) which satisfy the following relations From the relations above and the inverse matrix in eq. (A6), we obtain the following rotation rate Substituting eqs. (14) and (A11) into eq. (A4), we obtain the dynamics in the rotating frame where we have definedB 0 (t) = B 0 (t) − Γ rel . Expanding the first term on the right hand side and using the relations in eq. (A10) yields dR(t) −1 dt R(t) = −ωR I . One can rewrite eq. (A12) as where the stochastic operators transform asF where ∆ = Ω dc − ω with Ω dc = µ B g F B dc . The rotating frame transformation introduces terms oscillating at 2ω. If we apply apply the rotating wave approximation (RWA) by neglecting those terms such that Under the RWA the dynamical matrix M(t) would be then time independent M(t) = M rot . Regarding the contribution from external fields, the external matrix B ext (t) transforms into the matrix M ext (t) explicitly as where we have defined Ω y (t) = Ω ext y cos(ωt) + Ω ext x sin(ωt).
Hence, the generator of the dynamics can be defined asM Therefore, the dynamics in the rotating frame is written as where the stochastic operator transforms asF (t) = RF (t) and the generator of the dynamics can be defined according to eq. (A16) and according to eq. (A18), where we have defined Ω y (t) = Ω ext y cos(ωt) + Ω ext x sin(ωt).
The generatorM(t) is one of the most important matrices throughout the whole paper, since it contains all the interactions of the spin with the external magnetic fields and the relaxation terms.
At this point we can show that this result describes the magnetic resonance for polarized atomic spins. In this particular case, the relaxation matrix is Γ x = Γ y = Γ 2 and Γ z = Γ 1 Additionally, we consider no external magnetic field i.e. M ext (t) = 0, an average pumping rate Γ p (t) = Γ 0 and a polarized input state F in = (0, 0, F 0 z ), which is constant in the rotating frame R(t) F in = F in . Hence, according to eq. (A22), the steady state solution for the mean value is where the generator of the dynamics is defining Γ i = Γ i + Γ 0 with i = 1, 2. This solution can be explicitly written in the Cartesian components as where we have defined Ω rf = Ω rf /2. This result describes the magnetic resonance of a polarized atomic sample first obtained by Bloch [28].
1. Floquet expansion for the first moment of the spin operator in the rotating frame According to eq. (A22) the dynamics for the mean value of of the spin operator are given by where the stochastic noise contribution is F (t) = 0. One of the interesting features of describing the dynamics in the rotating frame, is that unlike the laboratory frame where the external fields enter as a constant variable, in the rotating frame the transverse external fields naturally exhibit a contribution to the first harmonics.
In particular, the external field can be decomposed as where with Ω ± = Ω ext x ± iΩ ext y . This result already shows that the transverse magnetic fields are mapped onto the quadrature of the spin evolution of the first harmonic e iωt since Ω ext x ± iΩ ext y . Similarly, the rotating frame matrix in eq. (A6) can be decomposed as where R (±1) = (R R ± iR I )/2 with the definition of R R and R I given in the Appendix A. This expansion implies that, according to the pump rate decomposition in eq. (19), the complete pumping term in eq. (A32), can be expanded as with In order to solve this harmonic dynamical equation, we employ the Floquet expansion for the spins in the rotating frame P (t) = n P (n) (t)e inωt . Applying the same procedure in Sec. IV B, we substitute the expansion P (t) and Γ (n) p (t) into eq. (A32) to obtain the the following dynamical equation for the vector of harmonics where with ext and the pump matrix term is (Γ in ) nm =Γ (n) p for n = m = 0, otherwise it is zero. The steady state solution takes the same form as in the laboratory frame case, therefore, However, the measurement of the spin dynamics is done in the laboratory frame. Therefore, applying the the inverse transformationF(t) = R(t)F (t), the steady state solution can be expressed in the harmonic linear space as where Appendix B: Diffusion matrix for spin operators The following description for Langevin dynamics is based in ref. [30]. The relaxation dynamics of the spin operators under the effect of stochastic operators can be written as We have defined in eq. (13) that the stochastic operatorŝ F i (t) describe a white noise process, which satisfy the following correlation function where F (t) = 0 andΓ ij = (Γ) ji . To determine the elements of the diffusion matrixΓ, we make use of the commutation relations of the spin operatorF i (t).
Let us start by solving the dynamical equation for the spin componentF i (t) from eqs. (B1-B3), such that Now, this solution must satisfy the commutation relation [F i (t),F j (t)] = i ijkFk (t) where ijk is the Levi-Civita tensor, for t and t + ∆t. Therefore, we can compute From the correlation function of the stochastic operators in eq. (B4) we obtain In particular, is worth noting that for ∆t = 0 the spin operators have no change, and therefore satisfy the commutator relations. Now, due to the decay process, for ∆t 1 the spin operator is slightly attenuated and compensated by the diffusion matrix. To notice that, let us examine the influence of the diffusion matrix into the short time dynamics, assuming ∆t 1, such that at first order If the diffusion matrix compensate the attenuation of the spin operator we would obtain the following solutioñ the spin operators would remain time independent, F k (t + ∆t) = F k (t) . The same procedure can be done for the anticommutator {F i (t),F j (t)}, such that In that particular case, from eqs. (B10) and (B11), one can find thatΓ where we have defined the second moment matrix in the lab frame as σ = F (t)F(t) T . The explicit expression for the diffusion matrix is However, in general, the diffusion matrix does not necessarily compensate the attenuation of the spin operators, it can reduce it though. Since the stochastic operators are modeling the flip in the atomic spins due to collision process, the second moment matrix can flip into the second moment of an unpolarized operator. Therefore, we can model the diffusion matrix as where σ 0 is the second moment matrix of an unpolarized sample that enters into the dynamics of the second moments of spin operators at the rate of Γ rel . Since σ 0 satisfy the commutation relations,Γ guarantees that the spin operators satisfy the commutation relations as is given in eq. (B8). The same procedure can be followed for the stochastic input operatorsF in (t). Considering that the three directions are equally pumped, the diffusion rate takes a simpler form from eq. (B14) where σ in correspond to the second moment matrix for an arbitrary input state.

Cross correlation functions
Let us consider a stochastic operatorŴ(t), where its mean value is Ŵ (t) = 0 and the correlation function is The time evolution of a symmetric correlation is From the correlation function in eq. (B16) we can write dt Γ w (t )δ(t − t). (B18) Using a change of variables, t = t − t (and t = t − t for the second integral) we obtain dt Γ w (t + t)δ(t ). (B19) In the case of a time independent diffusion matrix Γ w (t) = Γ w0 we find which is in agreement with ref. [32]. However for a time dependent diffusion matrix it is convenient to apply an > 0 around zero, for a proper definition of the integral with a δ(t) function, such that D w =2Γ w (t). (B21) The cross correlation of the spin operator with stochastic operators in eq. (32) is given in a general form as which generates a drift in the second moment dynamics proportional to the diffusion matrix Γ w (t). According to ref. [32], this matrix is non-zero becauseF(t) depends on the stochastic operators itself, therefore, the correlation W(t)F(t ) T = 0 when t = t, otherwise, there is no correlation since W(t) has a very short coherence time.
To determine the cross correlations D with either stochastic operators, W(t) ∈ {F in (t),F (t)}, let us consider from eq. (14), the time evolution of the spin operator as followŝ In the case where W(t) =F in (t), by multiplying from the right withF (t) T and taking the average we have It is worth pointing out that regardless the spin polarization F in , the third term is zero, since by definition F in (t) T = 0. Moreover, these kind of stochastic operators are only correlated with themselves, such that the correlation with any other operatorÔ(t) is Ô (t )F in (t) T = Ô (t ) F in (t) T , which by definition is zero. Therefore, the only non-zero term is Following the same procedure for F in (t)F ( t) T , and considering the case in eq. (B21) for the diffusion matrix in eq. (B15), the cross correlation with the input operator is Similarly for the unpolarized stochastic operatorF (t), the cross correlation is where the only non-zero term is Since the drift termΓ in eq. (13) is constant, according to eq. (B20), the complete cross correlation is From eq. (B14), we find that the drift matrix of the stochastic operator can take the form where σ 0 represents the second moment matrix of an unpolarized or thermal spin state.
Appendix C: Rotating Frame transformation of the second moment matrix dynamics Taking the time derivative of eq. (63) we have and from eq. (A11) we obtain Substituting eq. (35) into the rotating frame dynamics, considering the rotating wave approximation, leads to dσ(t) dt =M(t)σ (t) + σ (t)M(t) T − 2Γ p (t)σ (t) + Γ rel σ 0 + σ 0 Γ rel + 2Γ p (t) σ in (t) in whichM(t) is given in eq. (A22), the unpolarized second moment matrix remains diagonal and the relaxation matrix is such that In addition, the input second moment matrix is now time dependent σ in (t) = R(t) −1 σ in R(t) σ in (t) = σ in (0) + σ in (1) e iωt + σ in (−1) e −iωt , where σ in (0) =R (0) σ in R (0) Considering a non-polarized input state F in = 0, we can rewrite eq. (C3) as In the Liouville space, the dynamics equations for the second moment is given by eq. (64). In that case, the transformation of the input second moment matrix can be written as X in (t) = R in (t)X in with where R In the case where this transformation matrix is modulate by the pump rate, Γ p (t)R in (t), the effective time dependence can be harmonically expanded as Γ p (t) = Γ p (t)R in (t) = n Γ p (n) e inωt where Therefore, the dynamics in the Liouville space is given by 1. Transformation of the second moment matrix to the laboratory frame The rotating frame transformation for σ is Therefore, in the Liouville space, we can rewrite the transformation as X(t) = R(t)X (t) where R = L(R(t))R(R −1 (t)). Hence in terms of harmonics we have with where we have defined R (n) L = L(R (n) ) and R (n) R = R(R (n) ).