Exceptional points for parameter estimation in open quantum systems: analysis of the Bloch equations

Felix Bloch [1] pioneered the dynamical description of open quantum systems. Originally Bloch's equations describe the relaxation and dephasing of a nuclear spin in a magnetic field. Soon it became apparent that the treatment can be extended to a generic two-level-system (TLS), such as the dynamics of laser driven atoms in the optical regime [2–4]. The open TLS has been used to model many different fields of physics. The TLS or a q-bit is at the foundation of quantum information [5–9]. In particle physics the TLS algebra has been employed in studies of possible deviations from quantum mechanics in the context of neutrino oscillations [10], as well as quantum entanglement [11–15], associated with electron/positron collisions and entangled systems due to EPR-Bell correlations [16].

The TLS is the base for setting the frequency standard for atomic clocks [17]. As a result accurate measurement of frequency is an important issue. Quantum-enhanced measurements based on interferometry have been suggested as means to beat the shot noise limit [18]. In these methods the decoherence rate is the limiting factor [19]. In some cases quantum error correction can increase the coherence time and the accuracy [20]. In the present study we want to suggest an opposite strategy. By employing the non-Hermitian character of the dynamics, the decoherence can be transformed from a bug to a feature.

The Bloch equation is the simplest example of a quantum master equation. Bloch rederived the equation from first principles, employing the assumption of weak coupling between the system and bath [21, 22]. These studies have paved the way for a general theory of quantum open systems. Davies [23] rigorously derived the weak coupling limit, resulting in a quantum master equation which leads to a completely positive dynamical semigroup [24]. Based on a mathematical construction, Lindblad and Gorini, Kossakowski and Sudarshan (L-GKS) obtained the general structure of the generator ${\mathcal{L}}$ of a completely positive dynamical semigroup [25, 26]. In the Heisenberg representation the L-GKS generator becomes [27, ch 3]:

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}\hat{{\bf{X}}}\;=\;\displaystyle \frac{\partial \hat{{\bf{X}}}}{\partial t}+{\rm{i}}\left[\hat{{\bf{H}}},\hat{{\bf{X}}}\right]+\displaystyle \sum _{k}\left({\hat{{\bf{V}}}}_{k}^{\dagger }\hat{{\bf{X}}}{\hat{{\bf{V}}}}_{k}-\displaystyle \frac{1}{2}{\left[{\hat{{\bf{V}}}}_{k}^{\dagger }{\hat{{\bf{V}}}}_{k},\hat{{\bf{X}}}\right]}_{+}\right),\end{eqnarray} \tag{ 1 }$$

where $\hat{{\bf{X}}}$ is an arbitrary operator. The Hamiltonian $\hat{{\bf{H}}}$ is Hermitian and operators ${\hat{{\bf{V}}}}_{k}$ are defined to operate in the Hilbert space of the system. The [ · , ·] ₊ denotes an anti commutator.

The set of operators $\{\hat{{\bf{X}}}\}$ supports a Hilbert space construction using the scalar product: $({\hat{{\bf{X}}}}_{1},{\hat{{\bf{X}}}}_{2})\equiv \mathrm{tr}\{{\hat{{{\bf{X}}}_{1}}}^{\dagger }\hat{{{\bf{X}}}_{2}}\}.$ A crucial simplification to equation (1) is obtained when a set of operator is closed to the generator ${\mathcal{L}}.$ Then we can rephrase the dynamics with a matrix-vector notation [28]:

$$\begin{eqnarray}&&\dot{\vec{Y}}=M\vec{Y}\end{eqnarray} \tag{ 2 }$$

where $\vec{Y}$ is the vector of basis operators and M is the representation of the generator ${\mathcal{L}}$ in this vector space. The eigenvalues of the matrix M reflect the non-Hermitian dynamics generated by ${\mathcal{L}}.$ In general they are complex with the steady state eigenvector having an eigenvalue of zero. The solution for this equation is:

$$\begin{eqnarray*}&&\vec{Y}(t)={{\rm{e}}}^{{Mt}}\vec{Y}(0).\end{eqnarray*}$$

When M is diagonalizable, we can write M = T Λ T⁻¹, for a non-singular matrix T and a diagonal matrix Λ, which has the eigenvalues { λ_i } on the diagonal. Then we have e^{M t} = T e^{Λ t}T⁻¹ , with the diagonal matrix e^{Λ t}, which has the exponential of the eigenvalues, ${{\rm{e}}}^{{\lambda }_{i}t},$ on its diagonal. The resulting dynamics of expectation values of operators, as well as other correlation functions, follows a sum of decaying oscillatory exponentials. The analytical form of such dynamics is:

$$\begin{eqnarray}&&\left\langle X(t)\right\rangle =\displaystyle \sum _{k}\;{d}_{k}\mathrm{exp}[-{\rm{i}}{\omega }_{k}t],\end{eqnarray} \tag{ 3 }$$

where −i ω_k, denoted as complex frequencies, are the eigenvalues of M, d_k are the associated amplitudes, and both ω_k and d_k can be complex. The real part of the complex frequency ω_k represents the oscillation rate, while the imaginary part, $\mathrm{Im}({\omega }_{k})\leqslant 0$ represents the decaying rate.

For special values of the system parameters the spectrum of the non-Hermitian matrix M is incomplete. This is due to the coalescence of several eigenvectors, referred to as a non-Hermitian degeneracy. The difference between Hermitian degeneracy and non-Hermitian degeneracy is essential: in the Hermitian degeneracy, several different orthogonal eigenvectors are associated with the same eigenvalue. In the case of non-Hermitian degeneracy several eigenvectors coalesce to a single eigenvector [29, ch 9]. As a result, the matrix M is not diagonalizable.

The exponential of a non-diagonalizable matrix M can be expressed using its Jordan normal form: M = T J T⁻¹. Here, J is a Jordan-blocks matrix which has (at least) one non-diagonal Jordan block; J_i = λ_iI + N , where I is the identity and N has ones on its first upper off-diagonal. The exponential of M is expressed as e^Mt = T e^JtT⁻¹, with the block-diagonal matrix e^Jt, which is composed from the exponential of the Jordan blocks ${{\rm{e}}}_{i}^{J}t$. For non-Hermitian degeneracy of an eigenvalue λ_i, the exponential of the block J_i will have the form: ${{\rm{e}}}^{{J}_{i}t}={{\rm{e}}}^{{\lambda }_{i}{It}+{Nt}}\;=\;$ ${{\rm{e}}}^{{\lambda }_{i}t}{{\rm{e}}}^{{Nt}}.$ The matrix N is nilpotent and therefore the Taylor series of e^Nt is finite, resulting in a polynomial in the matrix Nt. This gives rise to a polynomial behaviour of the solution, and the dynamics of expectation values will have the analytical form of

$$\begin{eqnarray}&&\left\langle X(t)\right\rangle =\displaystyle \sum _{k}\displaystyle \sum _{\alpha =0}^{{r}_{k}}\;{d}_{k,\alpha }{t}^{\alpha }\mathrm{exp}[-{\rm{i}}{\omega }_{k}^{({r}_{k})}t],\end{eqnarray} \tag{ 4 }$$

replacing the form of equation (3). Here, ${\omega }_{k}^{({r}_{k})}$ denotes an eigenvalue with multiplicity of r_k + 1. Note that for non-degenerate eigenvalues, i.e. r_k = 0, we have d_k,0 = d_k and ω_k⁰ = ω_k. The difference in the analytic behaviour of the dynamics results in non-Lorentzian line shapes, with higher order poles in the complex spectral domain.

The point in the spectrum where the eigenvectors coalesce is known as an exceptional point (EP). When two eigenvalues of the master equation coalesce into one, a second-order non-Hermitian degeneracy is obtained. We refer to it as EP2, while a third-order non-Hermitian degeneracy is denoted by EP3.

This study addresses the scenario of the dynamics of a system coupled to a bath. The formalism is a reduced description of a tensor product of the system and the bath [27, 30]. The coupling to the bath introduces dissipation and dephasing into the dynamics. The state is represented as a density operator in Liouville space, and the dynamics is governed by the L-GKS equation. The non Hermitian properties of the dynamical generator ${\mathcal{L}}$ is caused by tracing out the bath degrees of freedom. We employ the Heisenberg picture with a complete operator basis set in Liouville space.

Previous studies of the physics of EPs investigated the scenario of scattering resonances phenomena. In that different scenario, the non Hermitian properties of the effective Hamiltonian are caused by the interaction between the discrete states via the common continuum of the scattering states [31, 32]. In those studies only coherent dynamics is considered and the dissipation and dephasing phenomena are absent.

Examples for EPs have been described in optics [33, 34], in atomic physics [35–40], in electron–molecule collisions [41], superconductors [42], quantum phase transitions in a system of interacting bosons [43], electric field oscillations in microwave cavities [44], in PT-symmetric waveguides [45], and in mesoscopic physics [46, 47].

Recently, Wiersig suggested a method to enhance the sensitivity of detectors using EPs [48]. Below we suggest to employ the EPs for the purpose of parameter estimation.

The analytical form of decaying exponentials, equation (3), is used in harmonic inversion methods to find the frequencies and amplitudes of the time series signal [49–51]. These frequencies and amplitudes can be employed to estimate the system parameters. If the sensitivity of the estimated frequencies is increased with respect to the system controls, the accuracy of the parameter estimation is enhanced. Such sensitivity increase can be achieved using the special character of the dynamics at EPs.

At EPs the analytical form includes also polynomials (equation (4)). Fuchs et al showed that applying the standard harmonic inversion methods to a signal generated by equation (4) leads to divergence of the amplitudes d_k. An extended harmonic inversion method can fix the problem. The divergence of the amplitudes d_k at the vicinity of EP can be used to locate them in the parameter space very accurately [52]. This is a consequence of the special non analytic character close to the EP (see in ch 9 in [29]).

Relying on the ability to accurately locate the EPs in the parameter space, we suggest to use the EPs for parameter estimation. The procedure we suggest follows:

• (i)  
  Accurately locate in the parameter space the desired EP by iterating the following steps:

  • (a)  
    Perform the experiment to get a time series of an observable for example the polarization as a function of time.
  • (b)  
    Obtain the characteristic frequencies and amplitudes of the signal using harmonic inversion methods.
  • (c)  
    In the parameter space, estimate the direction and distance to the EP and determine new parameters for the next iteration.
• (ii)  
  Invert the relations between the characteristic frequencies and the system parameters at the EP to obtain the system parameters.

The accurate location of the EPs, followed by inverting the relations, will lead to accurate parameter estimation.

4.1. The Bloch equation

The Bloch equation describes the dynamics of the three components of the nuclear spin, S_x, S_y, and S_z, under the influence of an external magnetic field, or a two-level atom in external electromagnetic field. In the rotating frame, we can write the equations in a matrix-vector notation:

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}\left(\begin{array}{c}{\tilde{S}}_{x}\\ {\tilde{S}}_{y}\\ {S}_{z}\end{array}\right)=\left(\begin{array}{ccc}-\displaystyle \frac{1}{{T}_{2}} & {\rm{\Delta }} & 0\\ -{\rm{\Delta }} & -\displaystyle \frac{1}{{T}_{2}} & \epsilon \\ 0 & -\epsilon & -\displaystyle \frac{1}{{T}_{1}}\end{array}\right)\left(\begin{array}{c}{\tilde{S}}_{x}\\ {\tilde{S}}_{y}\\ {S}_{z}\end{array}\right)+\left(\begin{array}{c}0\\ 0\\ \displaystyle \frac{1}{{T}_{1}}{S}_{z}^{0}\end{array}\right),\end{eqnarray} \tag{ 5 }$$

with T₁ and T₂ as the dissipation and dephasing relaxation parameters, and the detuning from resonance Δ and the amplitude as the field parameters. See details in appendix A.

The Bloch equations can be derived from the L-GKS equation of the two-level system, with the effective rotating-frame Hamiltonian

$$\begin{eqnarray*}&&\hat{{\bf{H}}}={\rm{\Delta }}{\hat{{\bf{S}}}}_{z}+\epsilon {\hat{\tilde{{\bf{S}}}}}_{x},\end{eqnarray*}$$

along with relaxation and dephasing terms. see appendix B for details.

Reducing the number of parameters, the master equation can be incorporated in the matrix:

$$\begin{eqnarray}{\rm{M}}=\left(\begin{array}{ccc}-\displaystyle \frac{{\rm{\Gamma }}}{2} & {\rm{\Delta }} & 0\\ -{\rm{\Delta }} & -\displaystyle \frac{{\rm{\Gamma }}}{2} & \epsilon \\ 0 & -\epsilon & -{\rm{\Gamma }}\end{array}\right),\end{eqnarray} \tag{ 6 }$$

with ${\rm{\Gamma }}=\frac{3}{2}\frac{1}{{T}_{1}}-\frac{1}{{T}_{2}}$ as the general relaxation coefficient (see appendix B).

The dynamics is determined by the exponential e^Mt, which typically describes oscillating decaying signal, see equation (3). Nevertheless, for specific parameters leading to EP the dynamics is modified to include polynomials, see equation (4).

4.2. EPs in the Bloch equation

The EPs are non-Hermitian degeneracies in the matrix M of equation (6). The task is to express the EPs using the parameters of this matrix. Explicit derivations are presented in appendix C. Non-Hermitian degeneracies of the eigenvalues [29], EP2, occur when

$$\begin{eqnarray*}&&{{\rm{\Gamma }}}^{4}{{\rm{\Delta }}}^{2}+16{\left({{\rm{\Delta }}}^{2}+{\epsilon }^{2}\right)}^{3}+{{\rm{\Gamma }}}^{2}\left(8{{\rm{\Delta }}}^{4}-20{{\rm{\Delta }}}^{2}{\epsilon }^{2}-{\epsilon }^{4}\right)=0.\end{eqnarray*}$$

Figure 1 shows a map of EP2 curve as a function of and Δ for fixed Γ = 0.1. Such figures were obtained in the study of analytical solutions for the Bloch equation [53–55].

Zoom In Zoom Out Reset image size

A third order EP, EP3, occurs when ${\rm{\Delta }}=\pm \sqrt{1/108}\;{\rm{\Gamma }},\epsilon =\sqrt{8/108}\;{\rm{\Gamma }}$ (red asterisks in figure 1). These triple-degeneracies EP3 occur twice, and have a cusp-like behaviour, emerging from the EP2-curves, identifiable as a section through an elliptic umbilic catastrophe [56]. This topology is also consistent with an analysis of non Hermitian degeneracies in a two-parameters family of 3 × 3 matrices [57]. In very strong driving fields the matrix M will loose symmetry [58, 59] maintaining the cusps but skewing the topology.

4.3. EP identification and parameter estimation

We now describe the two steps of the method for accurate determination the physical parameters. The first step is to identify the desired EP using a sequence of measured time-dependent signals. The second step is to invert the relations and determine the system parameters.

4.3.1. Identifying the second and third order EPs

To identify the EPs we used time series of the polarization observable ${S}_{z}\equiv \langle {\hat{{\bf{S}}}}_{z}\rangle ,$ initially at the ground state. We simulated the dynamics with varying field parameters (, Δ) generating a time series of polarization S_z[n] = S_z (n δ t). This signals served as the input for the harmonic inversion.

The parameters Δ and were tuned close to an EP. Generically we should have

$$\begin{eqnarray*}&&{S}_{z}(t)={d}_{1}{{\rm{e}}}^{-{\rm{i}}{\omega }_{1}t}+{d}_{2}{{\rm{e}}}^{-{\rm{i}}{\omega }_{2}t}+{d}_{3}{{\rm{e}}}^{-{\rm{i}}{\omega }_{3}t},\end{eqnarray*}$$

but in the EP2 (r_k = 1) we get

$$\begin{eqnarray*}&&{S}_{z}(t)={d}_{1}{{\rm{e}}}^{-{\rm{i}}{\omega }_{1}t}+\left({d}_{\mathrm{2,0}}+{d}_{\mathrm{2,1}}t\right){{\rm{e}}}^{-{\rm{i}}{\omega }_{2}^{(1)}t},\end{eqnarray*}$$

and for EP3 (r_k = 2)

$$\begin{eqnarray*}&&{S}_{z}(t)=\left({d}_{\mathrm{1,0}}+{d}_{\mathrm{1,1}}t+{d}_{\mathrm{1,2}}{t}^{2}\right){{\rm{e}}}^{-{\rm{i}}{\omega }_{1}^{(2)}t}.\end{eqnarray*}$$

(See equations (3) and (4).) We located suspected EPs by identifying possible degeneracies of the assigned frequencies ω_k. As stated earlier, applying standard harmonic inversion methods for the time series generated by a non-diagonalizable matrix, leads to divergence of the amplitudes d_k  [52]. This divergence can be used to locate the EPs accurately. A verification can be obtained by using the extended harmonic inversion method.

This procedure was employed to identify an EP2 for fixed Γ = 0.1 and  = 0.01, with varying Δ. The purple asterisks at figure 2 displays the absolute value of the difference between the frequencies $| {\omega }_{2}-{\omega }_{1}| ,$ obtained by the harmonic inversion for each parameter set. The degeneracy point is clearly observed. The diverging behaviour of the amplitudes is shown in red stars. It is consistent with the degeneracy of the frequencies. The EP2 is located at Δ = 1.021 × 10⁻³, consistent with the prediction. Using a finer mesh of sampling points the EP can be identified with a resolution exceeding 0.5 × 10⁻⁹.

Zoom In Zoom Out Reset image size

The EP3 was identified by a 2D search performed by varying and Δ, for fixed Γ = 0.1. We searched for the degeneracies of the three eigenvalues by employing the 2D function

$$\begin{eqnarray}&&F({\rm{\Delta }},\epsilon ,{\rm{\Gamma }})=\mathrm{log}\left(\left|\displaystyle \frac{1}{({\omega }_{1}-{\omega }_{2})}\displaystyle \frac{1}{({\omega }_{2}-{\omega }_{3})}\displaystyle \frac{1}{({\omega }_{3}-{\omega }_{1})}\right|\right),\end{eqnarray} \tag{ 7 }$$

which should diverges at the EP curve. Numerically, we get high values at this curve, with highest values obtained at the EP3. The upper panel of figure 3 shows the sharp curve of peaks following the curve of EPs. The highest point on the merging two ridges is the EP3. The lower panel of figure 3 shows the sum of the absolute values of the amplitudes, calculated by the standard harmonic inversion. The curve of the EPs is clearly identified.

Zoom In Zoom Out Reset image size

Refining the search leads to very high resolution, and the EP3 can be identified with a high accuracy, approaching the theoretical values of ${\rm{\Delta }}=\sqrt{1/108}\;{\rm{\Gamma }},$ $\epsilon =\sqrt{8/108}\;{\rm{\Gamma }}.$

An efficient algorithm to identify the EP3 is demonstrated based on a two-dimensional search in the parameter space of Δ and . This procedure enables the experimentalists to identify accurately the laser parameters for which the EP3 is obtained. We use the maximum of the function equation (7) as the objective leading to EP3.

Evaluating the function at each desired point in the parameter space include the following steps:

• (i)  
  Time series: obtain a time series of the polarization by performing the experiment or the numerical simulation.
• (ii)  
  Frequencies: calculate the frequencies from the time series by harmonic inversion.
• (iii)  
  Function evaluation: evaluate the function F(Δ, , Γ) from the calculated frequencies.

Standard search methods can stagnate due to the high values at the EP2 curve. Another difficulty is the cusp behaviour of the EP2 curve close to the EP3. To overcome these diffculties we implemented a 'climbing the valley' procedure: staying on the valley of the local minima ensures the search overcomes the stagnation due to the EP2 curve. The procedure follows:

• (i)  
  Preliminary step—initial point:

  • (a)  
    Locate points inside the triangle-like EP curve (see figure 4). The inner area of the curve is characterized by real-only eigenvalues.
  • (b)  
    Perform a 1D search to find a minimum on a straight line.
• (ii)  
  Valley ascend: each iteration ascends up the valley to a point with higher value of the function F. This is done by finding a minimum on the circular arc that is centred at the current point, enclosed by two radii. The angles of these radii can be predefined or defined on each iteration. We perform the following steps:

  • (a)  
    Determining the angular range. Predefined or from the previous iterations.
  • (b)  
    Determining the radius. The radius is the distance from the current point to nearest point on the EP2 curve that is in the angular range.
  • (c)  
    Finding the next point. Performing a 1D search on the circular arc that is defined by the angular range and the radius (see blue arc in figure 4). The point for the next iteration is the point on the arc with the minimal value of F (see end of green line in figure 4).

Zoom In Zoom Out Reset image size

These steps converge to the desired EP3 point. Figure 4 demonstrates the progress in the 'valley ascend' method with a few iterations.

The Valley ascend method presented above is a generic method, and can be used also for searching higher order degeneracies in other systems. For the Bloch equation case, where the generating matrix, equation (6), is a 3 × 3 matrix, the EP3 is the point where the characteristic polynomial

$$\begin{eqnarray}&&{P}_{{}_{{\rm{\Delta }},\epsilon ,{\rm{\Gamma }}}}(\omega )=(\omega -{\omega }_{1})(\omega -{\omega }_{2})(\omega -{\omega }_{3})\end{eqnarray} \tag{ 8 }$$

has roots with multiplicity of 3. Therefore we can use the special properties of the cubic equation and perform a regular root search. We define r, s and t as the coefficient of the polynomial P_{Δ, , Γ}(ω) defined in equation (8):

$$\begin{eqnarray}&&(\omega -{\omega }_{1})(\omega -{\omega }_{2})(\omega -{\omega }_{3})={\omega }^{3}+r{\omega }^{2}+s\omega +t.\end{eqnarray} \tag{ 9 }$$

We define the functions

$$\begin{eqnarray}p({\rm{\Delta }},\epsilon ,{\rm{\Gamma }}) & = & s-\displaystyle \frac{1}{3}{r}^{2}\\ q({\rm{\Delta }},\epsilon ,{\rm{\Gamma }}) & = & \displaystyle \frac{2}{27}{r}^{3}-\displaystyle \frac{1}{3}{rs}+t,\end{eqnarray} \tag{ 10 }$$

and perform a 2D conventional root search. The point in the parameter space where these two functions vanish is point where the three eigenvalues are degenerate. We have applied this method using standard method of 2D root search obtaining high accurate values of the EP3.

4.3.2. Physical parameter estimation from the value at the EP

For the TLS the system parameters are the frequency ω_s associated with the energy gap, the general decay rate Γ and the dipole strength μ. The external experimentally controlled parameters are the driving frequency ν and the power amplitude ${\mathcal{E}}$. The parameters of equation (6) can be related with $\epsilon =\mu {\mathcal{E}}$ and Δ = ω_s − ν. One would like to estimate the system parameters from experiments. After locating accurately the EP, we can determine the parameters by inverting the relations between the eigenvalues and the system parameters.

To obtain high accuracy, we used the identification the triple-degeneracy point EP3 presented above, so both parameters—Δ and — are located accurately. The accurate location of Δ and makes the parameter estimation very robust to uncertainties in the location of the EPs. This is a consequence of the special non analytic character close to the EP3 (see appendix D). Therefore, the system parameters Γ, ω_s and μ can be determined to a high degree of accuracy at this point. From the eigenvalues of the matrix M in equation (6) we get ${\rm{\Gamma }}=\frac{{\rm{i}}}{2}\left({\omega }_{1}+{\omega }_{2}+{\omega }_{3}\right).$ To obtain and Δ one has to invert nonlinear relations (see appendix C). At the EP3, the inversion becomes: ${\omega }_{s}=\nu +\sqrt{1/108}\;{\rm{\Gamma }},$ $\mu =\sqrt{8/108}\;{\rm{\Gamma }}/{\mathcal{E}}.$

4.3.3. Noise sensitivity

Parameters estimation naturally raises the issue of sensitivity to noisy experimental data. The noise sensitivity will be determined by the method of harmonic inversion. If the sampling periods have high accuracy then the time series can be shown to have an underlying Hamiltonian generator. This is the basis for linear methods, such as the filter diagonalization (FD) [49, 50]. The noise in these methods results in normally distributed underlying matrices, and the model displays monotonous behaviour with respect to the noise. This was verified analytically and by means of simulations in [60]. As a result sufficient averaging will eliminate the noise. Practical implementations require further analysis with evidence of nonlinear effects of noise. For example, Mandelshtam et al analysed the noise-sensitivity of the FD in the context of NMR experiments [61, 62] and Fourier transform mass spectrometry [63]. For some other methods, a noise reduction technique was proposed in [51].

Bloch's equation has become the template for the dynamics of open quantum systems. Such systems typically decohere with a dynamical signature of decaying oscillatory motion. It is therefore surprising that the existence of non Hermitian degeneracies has been overlooked. Our finding of an intricate manifold of double degeneracies EP2 and triple degeneracies EP3 in the elementary TLS template suggests that any quantum dynamics described by the L-GKS generator [25, 26] will exhibit a manifold of EPs.

Non Hermitian degeneracies of the EP have a subtle influence on the dynamics. The hallmark of EP dynamics is a polynomial component in the decay leading to non-Lorentzian lineshapes. We suggest an experimental procedure to identify the EP in Bloch systems, using harmonic inversion of the polarization time series. The sensitivity of harmonic inversion in the neighbourhood of an EP enables us to accurately locate the EP, and therefore allows us to determine the system parameters: the energy gap ω_s, the dipole transition moment μ, and the decoherence rate Γ.

This study is only the first step in establishing parameter estimation via EPs. A generalization to larger Liouville spaces is under study for atomic spectroscopy. Under the influence of driving fields and due to spontaneous emission, atoms and ions can have a structure of N-level system with relaxation. In these systems we expect non-Hermitian degeneracy of high order. The structure of the EPs in these systems can be used for estimating the energy differences, the lifetimes, and branching ratios. Work in this direction is in progress.

Many quantum systems are open and their dynamics has dissipative nature, which is described well by the L-GKS equation. Therefore we expect to find EPs in many quantum systems. Under the appropriate circumstances these EPs can be used for accurate parameter estimation.

We thank Ido Schaefer, Amikam Levy, and Raam Uzdin for fruitful discussions. We thank Jacob Fuchs and Jörg Main for assisting with the extended harmonic inversion method. We thank the referee for proposing the root search for the EP3. Work supported by the Israel Science Foundation Grants No. 2244/14 and No. 298/11 and by I-Core: the Israeli Excellence Center 'Circle of Light'.

The Bloch equation describes the dynamics of the three components of the nuclear spin, S_x, S_y, and S_z, under the influence of an external magnetic field $\vec{H}.$ The equations as appear in Bloch's original paper ([1], equation (38)) are

$$\begin{eqnarray}\dot{{S}_{x}} & = & \gamma \left({S}_{y}{H}_{z}-{S}_{z}{H}_{y}\right)-\displaystyle \frac{1}{{T}_{2}}{S}_{x}\\ \dot{{S}_{y}} & = & \gamma \left({S}_{z}{H}_{x}-{S}_{x}{H}_{z}\right)-\displaystyle \frac{1}{{T}_{2}}{S}_{y}\\ \dot{{S}_{z}} & = & \gamma \left({S}_{x}{H}_{y}-{S}_{y}{H}_{x}\right)-\displaystyle \frac{1}{{T}_{1}}\left({S}_{z}-{S}_{z}^{0}\right).\end{eqnarray} \tag{ A.1 }$$

T₁ and T₂ are two relaxation parameters (the pure dephasing rate $\frac{1}{{T}_{2}^{*}}$ is related by $\frac{1}{{T}_{2}}=\frac{1}{2{T}_{1}}+\frac{1}{{T}_{2}^{*}}$), γ is the gyromagnetic ratio, and S_z⁰ is the equilibrium value of S_z under the influence of constant external magnetic field H_z = H₀. These equations can be recast in a matrix-vector notation:

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}\left(\begin{array}{c}{S}_{x}\\ {S}_{y}\\ {S}_{z}\end{array}\right)=\left(\begin{array}{ccc}-\displaystyle \frac{1}{{T}_{2}} & \gamma {H}_{z} & -\gamma {H}_{y}\\ -\gamma {H}_{z} & -\displaystyle \frac{1}{{T}_{2}} & \gamma {H}_{x}\\ \gamma {H}_{y} & -\gamma {H}_{x} & -\displaystyle \frac{1}{{T}_{1}}\end{array}\right)\left(\begin{array}{c}{S}_{x}\\ {S}_{y}\\ {S}_{z}\end{array}\right)+\left(\begin{array}{c}0\\ 0\\ \displaystyle \frac{1}{{T}_{1}}{S}_{z}^{0}\end{array}\right).\end{eqnarray} \tag{ A.2 }$$

For an external field $\vec{H}$ with the components ${H}_{x}={H}_{1}\mathrm{cos}\omega t,$ ${H}_{y}=-{H}_{1}\mathrm{sin}\omega t$, H_z = H₀, we define the rotating frame:

$$\begin{eqnarray}{S}_{x} & = & \tilde{{S}_{x}}\mathrm{cos}\omega t-\tilde{{S}_{y}}\mathrm{sin}\omega t\\ {S}_{y} & = & -\tilde{{S}_{x}}\mathrm{sin}\omega t-\tilde{{S}_{y}}\mathrm{cos}\omega t.\end{eqnarray} \tag{ A.3 }$$

With the notations  = γ H₁ and Δ = γ H₀ − ω we have (see also [4]):

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}\left(\begin{array}{c}{\tilde{S}}_{x}\\ {\tilde{S}}_{y}\\ {S}_{z}\end{array}\right)=\left(\begin{array}{ccc}-\displaystyle \frac{1}{{T}_{2}} & {\rm{\Delta }} & 0\\ -{\rm{\Delta }} & -\displaystyle \frac{1}{{T}_{2}} & \epsilon \\ 0 & -\epsilon & -\displaystyle \frac{1}{{T}_{1}}\end{array}\right)\left(\begin{array}{c}{\tilde{S}}_{x}\\ {\tilde{S}}_{y}\\ {S}_{z}\end{array}\right)+\left(\begin{array}{c}0\\ 0\\ \displaystyle \frac{1}{{T}_{1}}{S}_{z}^{0}\end{array}\right).\end{eqnarray} \tag{ A.4 }$$

These equations also describe, in the dipole approximation, a two-level atom in external electromagnetic field. In this case, the system parameters are the the unperturbed frequency of the system ω_s, and the dipole strength μ. The external experimentally controlled parameters are the driving frequency ν and the power amplitude ${\mathcal{E}}.$ The parameters of equation (A.4) are related with $\epsilon =\mu {\mathcal{E}}$ and Δ = ω_s − ν. In the absence of dissipation the eigenvalues of the matrix are pure imaginary, and the dynamics is a free precession of the polarization vector characterized by the Rabi frequency: ${\rm{\Omega }}=\sqrt{{\epsilon }^{2}+{{\rm{\Delta }}}^{2}}.$ When dissipation is present the eigenvalues of the homogeneous part of equation (A.4) become complex, reflecting a decaying oscillation dynamics leading asymptotically to a steady state.

In the Heisenberg representation the L-GKS generator becomes:

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}\hat{{\bf{X}}}\;=\;\displaystyle \frac{\partial \hat{{\bf{X}}}}{\partial t}+{\rm{i}}\left[\hat{{\rm{H}}},\hat{{\bf{X}}}\right]+\displaystyle \sum _{k}\left({\hat{{\rm{V}}}}_{k}^{\dagger }\hat{{\bf{X}}}{\hat{{\rm{V}}}}_{k}-\displaystyle \frac{1}{2}\left\{{\hat{{\rm{V}}}}_{k}^{\dagger }{\hat{{\rm{V}}}}_{k},\hat{{\bf{X}}}\right\}\right).\end{eqnarray} \tag{ B.1 }$$

where $\hat{{\bf{X}}}$ is an arbitrary operator. The Hamiltonian $\hat{{\bf{H}}}$ is Hermitian and $\hat{{\bf{V}}}$ is defined to operate in the Hilbert space of the system. The curly brackets denote an anti commutator. The set of operators $\{\hat{{\bf{X}}}\}$ supports a Hilbert space construction, with the scalar product defined as: $\left({\hat{{\bf{X}}}}_{1},{\hat{{\bf{X}}}}_{2}\right)\equiv \mathrm{tr}\left\{{\hat{{\bf{X}}}}_{1}^{\dagger }{\hat{{\bf{X}}}}_{2}\right\}.$

For two-level system, the effective rotating-frame Hamiltonian under a driving field with detuning Δ and driving frequency is:

$$\begin{eqnarray}&&\hat{{\bf{H}}}={\rm{\Delta }}{\hat{{\bf{S}}}}_{z}+\epsilon {\tilde{{\rm{S}}}}_{x}.\end{eqnarray} \tag{ B.2 }$$

The TLS L-GKS equation for an operator $\hat{{\bf{X}}}$ with relaxation and pure dephasing becomes

$$\begin{eqnarray}\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}\hat{{\bf{X}}} & = & {\rm{i}}\;[\hat{{\bf{H}}},\hat{{\bf{X}}}]\\ & & +{\kappa }_{-}\left({\hat{{\bf{S}}}}_{+}\hat{{\bf{X}}}{\hat{{\bf{S}}}}_{-}-\displaystyle \frac{1}{2}\{{\hat{{\bf{S}}}}_{+}{\hat{{\bf{S}}}}_{-},\hat{{\bf{X}}}\}\right)\\ & & +{\kappa }_{+}\left({\hat{{\bf{S}}}}_{-}\hat{{\bf{X}}}{\hat{{\bf{S}}}}_{+}-\displaystyle \frac{1}{2}\{{\hat{{\bf{S}}}}_{-}{\hat{{\bf{S}}}}_{+},\hat{{\bf{X}}}\}\right)\\ & & -\gamma \;[{\hat{{\bf{S}}}}_{z},[{\hat{{\bf{S}}}}_{z},\hat{{\bf{X}}}]],\end{eqnarray} \tag{ B.3 }$$

where κ_± are kinetic coefficients, ${\kappa }_{+}/{\kappa }_{-}=\mathrm{exp}(-{\hslash }\omega /{k}_{B}T),$ and γ is the pure dephasing rate [2, 64].

To rephrase the equation in a matrix-vector notation, We use the polarization operators and the identity matrix to form the vector of basis operators: ${\vec{S}}^{\prime }={\left({\tilde{{\bf{S}}}}_{x},{\tilde{{\bf{S}}}}_{y},{\hat{{\bf{S}}}}_{z},\hat{{\bf{I}}}\right)}^{T}.$ Then equation (B.3) can be written as $\dot{{\vec{S}}^{\prime }}={M}^{\prime }{\vec{S}}^{\prime },$ with an appropriate 4 × 4 matrix M^'. We can reduce the dimensions by writing an inhomogeneous equation for the three-component vector $\vec{S}={\left({\tilde{{\bf{S}}}}_{x},{\tilde{{\bf{S}}}}_{y},{\hat{{\bf{S}}}}_{z}\right)}^{T}:$

$$\begin{eqnarray}&&\dot{\vec{S}}=\left(M-\gamma I\right)\left(\vec{S}-{\vec{S}}_{\mathrm{eq}}\right),\end{eqnarray} \tag{ B.4 }$$

with Γ = κ₋ + κ₊ − γ, I as the 3 × 3 identity matrix, ${\vec{S}}_{\mathrm{eq}}$ that fulfills $(\gamma I-M){\vec{S}}_{\mathrm{eq}}={(0,\;0,({\kappa }_{+}-{\kappa }_{-})\hat{{\bf{I}}})}^{T}$ and the matrix:

$$\begin{eqnarray}{\rm{M}}=\left(\begin{array}{ccc}-\displaystyle \frac{{\rm{\Gamma }}}{2} & {\rm{\Delta }} & 0\\ -{\rm{\Delta }} & -\displaystyle \frac{{\rm{\Gamma }}}{2} & \epsilon \\ 0 & -\epsilon & -{\rm{\Gamma }}\end{array}\right).\end{eqnarray} \tag{ B.5 }$$

Equation (B.4) can be merged with the Bloch's equation (A.4) where $\frac{1}{{T}_{1}}={\kappa }_{+}+{\kappa }_{-}$ and $\frac{1}{{T}_{2}}\;=\;$ $\gamma +\frac{1}{2}({\kappa }_{+}+{\kappa }_{-}).$

The general solution for this equation is:

$$\begin{eqnarray}&&\vec{S}(t)={{\rm{e}}}^{-\gamma t}{{\rm{e}}}^{{Mt}}({\vec{S}}_{0}-{\vec{S}}_{\mathrm{eq}})+{\vec{S}}_{\mathrm{eq}},\end{eqnarray} \tag{ B.6 }$$

with ${\vec{S}}_{0}=\vec{S}(0).$

The master equation equation (B.3) is a common form for TLS found in the literature [7, 65, 66]. Equation (B.5) which determines the EP interpolates between two extreme cases. The first is associated with spontaneous emission, then Γ = κ₋ . The second is a hot singular bath dominated by pure dephasing, then Γ =  − γ.

The task is to find the eigenvalues of the generator matrix (6).

We first define the variables:

$$\begin{eqnarray}Y & = & 12{{\rm{\Delta }}}^{2}+12{\epsilon }^{2}-\ {{\rm{\Gamma }}}^{2}\\ X & = & -36{{\rm{\Delta }}}^{2}+18{\epsilon }^{2}-{{\rm{\Gamma }}}^{2}.\end{eqnarray} \tag{ C.1 }$$

We also define:

$$\begin{eqnarray}W & = & \sqrt{{{\rm{\Gamma }}}^{2}{X}^{2}+{Y}^{3}}\\ & = & {\left({{\rm{\Gamma }}}^{4}{{\rm{\Delta }}}^{2}+16{\left({{\rm{\Delta }}}^{2}+{\epsilon }^{2}\right)}^{3}+{{\rm{\Gamma }}}^{2}\left(8{{\rm{\Delta }}}^{4}-20{{\rm{\Delta }}}^{2}{\epsilon }^{2}-{\epsilon }^{4}\right)\right)}^{1/2}.\end{eqnarray} \tag{ C.2 }$$

With these definitions the eigenvalues of equation (6) become:

$$\begin{eqnarray}{m}_{1} & = & -\displaystyle \frac{2}{3}{\rm{\Gamma }}+\displaystyle \frac{1}{6}\left({\left(W+{\rm{\Gamma }}X\right)}^{1/3}-\displaystyle \frac{Y}{{\left(W+{\rm{\Gamma }}X\right)}^{1/3}}\right)\\ {m}_{2} & = & -\displaystyle \frac{2}{3}{\rm{\Gamma }}+\displaystyle \frac{1}{6}\left({{\rm{e}}}^{{\rm{i}}\frac{2}{3}\pi }{\left(W+{\rm{\Gamma }}X\right)}^{1/3}+{{\rm{e}}}^{{\rm{i}}\frac{1}{3}\pi }\displaystyle \frac{Y}{{\left(W+{\rm{\Gamma }}X\right)}^{1/3}}\right)\\ {m}_{3} & = & -\displaystyle \frac{2}{3}{\rm{\Gamma }}+\displaystyle \frac{1}{6}\left({{\rm{e}}}^{-{\rm{i}}\frac{2}{3}\pi }{\left(W+{\rm{\Gamma }}X\right)}^{1/3}+{{\rm{e}}}^{-{\rm{i}}\frac{1}{3}\pi }\displaystyle \frac{Y}{{\left(W+{\rm{\Gamma }}X\right)}^{1/3}}\right).\end{eqnarray} \tag{ C.3 }$$

For real W (i.e. for ${{\rm{\Gamma }}}^{2}{X}^{2}+{Y}^{3}\geqslant 0$) all eigenvalues are real. For Γ²X² + Y³ < 0, W is complex, and two of the eigenvalues are complex (complex conjugate to each other).

Non-Hermitian degeneracies of the eigenvalues occur when W vanishes. In such cases the second and third eigenvalues are degenerated, leading to EP2. A third order EP, EP3, occurs for X = Y = 0. This happens when ${\rm{\Delta }}=\pm \sqrt{1/108}\;{\rm{\Gamma }},\epsilon =\sqrt{8/108}\;{\rm{\Gamma }}.$ These triple-degeneracies EP3 occur twice, and have a cusp-like behaviour, emerging from the EP2-curves, identifiable as an elliptic umbilic catastrophe [56]. This topology is also consistent with the analysis of non Hermitian degeneracies of a two-parameters family of 3 × 3 matrices, done by Mailybaev [57]. In very strong driving fields the matrix M will loose symmetry [58, 59] maintaining the cusps but skewing the topology.

There is a special non analytic character close to the EP3: when $\nu \to {\nu }^{\mathrm{EP}3}$ and ${\mathcal{E}}\to {{\mathcal{E}}}^{\mathrm{EP}3}$ then the three frequencies obtained by the standard harmonic inversion coalesce, leading to a branch point (see ch 9 in [29]):

$$\begin{eqnarray}&&{\omega }_{k=\mathrm{1,2,3}}={\omega }_{1}^{(2)}+{{\rm{e}}}^{{\rm{i}}\frac{2\pi }{3}}{\left[{\alpha }_{k}(\nu -{\nu }^{\mathrm{EP}3})+{\beta }_{k}({\mathcal{E}}-{{\mathcal{E}}}^{\mathrm{EP}3})\right]}^{\frac{1}{3}}\end{eqnarray} \tag{ D.1 }$$

where α_k and β_k are parameters. At the EP3, i.e. for $\nu \to {\nu }^{\mathrm{EP}3}$ and ${\mathcal{E}}\to {{\mathcal{E}}}^{\mathrm{EP}3},$ we get $\partial {\omega }_{k}/\partial \nu \to \infty$ and $\partial {\omega }_{k}/\partial {\mathcal{E}}\to \infty ,$ leading to $\partial {\rm{\Gamma }}/\partial \nu \to \infty$ and $\partial {\rm{\Gamma }}/\partial {\mathcal{E}}\to \infty .$