Quantum Fisher information of two atoms with dipole–dipole interaction under the environment of phase noise lasers

We investigate the parameter estimation problems of two-atom system driven by the phase noise lasers (PNLs) environment. And we give a general method of numeric solution to handle the problems of atom system under the PNLs environment. The calculation results of this method on Quantum Fisher Information (QFI) are consistent with our former results. Moreover, we consider the dipole–dipole (d–d) interaction between the atoms under PNLs environment with the collective decay, and the results show that larger d–d interaction and smaller collective decay rate lead to larger QFI of the two-atom system. So the collective decay will destroy the QFI while the d–d interaction will preserve the QFI, these results can be used to protect the QFI of two-atom system driven by the PNLs environment.

This equation is the general expression of QFI in the diagonalized space of the density matrix ρ(φ). Now we begin to discuss the situation where the atoms have the qubit-qubit interactions. For simplicity, we only discuss two qubits (or atoms) interact with two individual PNLs 20,21 , i.e., each of the two atoms is assumed to be driven by a classical field with a randomly fluctuating phase � l (t) , where l = 1, 2 indicates the two atoms separately. They are the two lasers' phase noises respectively which act on their respective atoms by the interaction terms. This noise � 1 (t) (or � 2 (t) ) means laser's phase will diffuse with time, and the initial phase will no longer be a constant. So the Hamiltonian with the interaction part of the atom and with the phase noise � l (t) of the PNL is where one of the two atoms with frequency ω 0l is described by the upper state |0� and ground state |1� with σ .
. www.nature.com/scientificreports/ σ 1(2) z = |0� 1(2) �0| − |1� 1(2) �1| = 1 0 0 − 1 ⊗ I 2 (I 2 ⊗ 1 0 0 − 1 ) 41 . As we know, if the phase noises are random telegraph noise (RTN), the RTN model means that � 1 (t) (or � 1 (t) ) jumps between two values, e.g., a and −a 1,44,45 , which is a Poisson process. Here we discuss laser phase noise of the standard phase diffusion model, the phase noise � l (t) is a simple one-dimensional random walk or Wiener process, i.e., � l (t) is white noise: ��� l (t)� l (t + τ )�� = 2D l δ(τ ) , l = 1, 2 . So the atom is effectually being driven by complex colored noise e i�(t) with correlation function where subscripts l = 1, 2 indicate the two atoms separately. It is the fact that the driving term is colored noise that ultimately leads to a non-Markovian master equation. Let us consider the Hamiltonian (9) in the rotating frame. Here we use the unitary transformation with operator U = exp (− i 2 2 l=1 ω lσ l z t) to act in the Hamiltonian (9). And under rotating-wave approximation, in the basis of 4 × 4 : E (2) = two qubits' Hamiltonian in the rotating frame (i.e., under the unitary transformation of U) is Here we have chosen the same for the two qubits, and the two qubits are assumed to have the d-d interaction (or qubit-qubit interaction) with the coefficient g. l = ω 0l −ω l (l = 1, 2) is the detuning and we assume the two atoms are driven by similar noise lasers (i.e., ω 1 = ω 2 = ω ). Below we let 1 = 2 = 0 for simplicity. So we have ω 01 = ω 02 = ω 0 = ω , which means the transition frequency of the two atoms is resonance with the lasers. For calculating the QFI, our purpose is to obtain the two qubits' last evolution state ρ(φ) which is the stochastic average 39 of all the possible evolution states. When the atoms are closer, there is collective decay 2,3,40 , and here we use the master equation with collective decay 39 ± . Note that the coefficient matrix (γ ij ) ij of the collective decay has not only diagonal but also non-diagonal elements, which are the spontaneous decay rates for the cooperative system 2,40 .
Let us consider how to construct the probability distribution P(� 1 , � 2 , t) . If we assume that in a small time interval t , there is equal probability R l t of l increasing or decreasing by an amount l then we have In the limit of t → 0 , is the constant diffusion rates of the two PNLs, respectively. And Eq. (14) leads to the well-known probability diffusion equation 31 Below we will derive the master equation of quantum density operator ρ(� 1 , � 2 , t) under the evolution of random phase noise l , l = 1, 2 . It will be normalized to unity, i.e., Tr[ρ(� 1 , � 2 , t)] = 1 . The density operator ρ(� 1 , � 2 , t) will be the density operator obtained by averaging over all the previous history of evolution that lead to the random variable l having the current value l .For small time interval t and small phase intervals: 1 and 2 , the density operator ρ(� 1 , � 2 , t + �t) will be made up from five contributions depending on whether the ρ(� 1 , � 2 , t + �t) evolves freely with no change in 1 and 2 with the probability of (1 − (R 1 + R 2 )�t) , or else there have occurred four 'jumps' 31,38 1 2 R 1 t for each one, and ρ(� 1 , � 2 ± �� 2 , t) → ρ(� 1 , � 2 , t + �t) with the probability of 1 2 R 2 t for each one. The states ρ(� 1 , � 2 , t) , ρ(� 1 ± �� 1 , � 2 , t) and ρ(� 1 , � 2 ± �� 2 , t) have to be separately weighted by the probabilities P(� 1 , � 2 , t) , P(� 1 ± �� 1 , � 2 , t) and P(� 1 , � 2 ± �� 2 , t) at time t. Putting this all together and normalizing the state yield where Here we have considered the zero-point fluctuations and 2 γ is equal to Einstein A coefficient for a single atom and S ± = 2 l=1 S l ± . In the limit of t → 0 , 1 → 0 , 2 → 0 , R 1 → ∞ and R 2 → ∞ , by using Eqs. (14) and (15), Eq. (17) can be transformed into We assume that the phase 1 and 2 is initially uniformly distributed over [0, 2π ), and hence for all time, the probability distribution of the phase 1 and 2 will be its steady state distribution P( � 1 , � 2 , t) = 1 2π . And we assume the initial condition ρ( � 1 , � 2 , 0) = ρ( 0) , which means the system is brought into interaction with the stochastic influence in a state that is uncorrelated with the stochastic driving term H(� 1 , � 2 ) . As the P( � 1 , � 2 , t) is independent of time, Eq. (19) will be simplified to If we discretize the ρ(� 1 , � 2 , t) on the stochastic phases 1 , 2 and time t with discrete values, we can give the relation of subscripts i 1 and i 2 with 1 and 2 : � 1 = (i 1 − 1)�� 1 and � 2 = (i 2 − 1)�� 2 . Then the matrices have three kinds of equal expressions by using subscripts i 1 and i 2 , 1 and 2 and (i 1 − 1) 1 and (i 2 − 1)�� 2 respectively: ρ n i 1 , . And these three kind expressions are the same here and after even under later discussions where these matrices are rearranged. So Eq. (20) can be expressed by the finite difference scheme (here we use the implicit scheme of the finite difference scheme): In this equation, n is the time discrete index and can be any value of 1, 2, . . . , n m ( n m is the maximal time index). And for the superscript, ρ n means the solution of the density matrix at n (or the time t), and ρ n m means the solution of ρ n at n = n m (or the final time t = t m ). And the density matrix ρ n−1 is the solution of the density matrix at n − 1 (or the any time t − t ), which means ρ n−1 is the one that t before ρ n . Note that ρ 0 means initial time value of the density matrix ρ n−1 at n = 1 (or the moment t = 0 ). And ρ 1 means initial time value of (17) www.nature.com/scientificreports/ the density matrix ρ n at n = 1 (which means the moment t = t ), and also it can be matrix ρ n−1 at n = 2 (which also means the moment t = t ). This is the base for our iteration calculation later. And �t = (t m − 0)/n m which is the time step length. Here �� 1 = (2π − 0)/N 1 and �� 2 = (2π − 0)/N 2 are the step lengths of diffusion of phase noise, and N 1 and N 2 are maximal phase index respectively (i.e., i 1 and i 2 have the maximal values N 1 and N 2 respectively). In Eq. (21), for the subscripts, we further choose i 1 and i 2 from 1 to N 1 and 1 to N 2 as the phase noise indexes, which means the possible values of phase noise 1 and 2 , respectively. And then we can get a series of equations, namely, N 1 N 2 equations. Summarizing these N 1 N 2 equations, we can get a matrix equation: Note that in the summarizing process we have used ρ n 0,i 2 = ρ n N 1 ,i 2 and ρ n i 1 ,0 = ρ n i 1 ,N 2 , because of ρ(−�� 1 , . Thus we can ensure subscripts satisfying 1 ≤ i 1 ≤ N 1 and 1 ≤ i 2 ≤ N 2 . Here is an unit matrix with its dimension of N l × N l here and after. I ⊗2 2 is the unit operator of two qubits to each block. This matrix equation canbe simply expressed as A n d h e re we h ave g ive n a s i mpl i f i c at i on :  , in which the two adjacent diagonals of the main diagonal is 1, the matrix element of the top right and bottom left is 1, and other elements are 0. Here X = (ρ n i 1 ,i 2 ) 1≤i 1 ≤N 1 ,1≤i 2 ≤N 2 and C = (ρ n−1 i 1 ,i 2 ) 1≤i 1 ≤N 1 ,1≤i 2 ≤N 2 (below we express ( ) 1≤i 1 ≤N 1 ,1≤i 2 ≤N 2 as ( ) i 1 ,i 2 for simplicity). X and C are both N 1 × N 2 block matrices displayed by subscript ( ) i 1 ,i 2 , and each block is 4 × 4 matrix. Here the first value on moment zero is given by ρ(� 1 , � 2 , 0) = (ρ 0 i 1 ,i 2 ) 1≤i 1 ≤N 1 ,1≤i 2 ≤N 2 , which means the first value of n is 0. And the last value on moment t m is given by ρ(� 1 , � 2 , t m ) = (ρ n m i 1 ,i 2 ) 1≤i 1 ≤N 1 ,1≤i 2 ≤N 2 , which means (n m + 1) th value of n is n m .
. Here X and C are the matrixes to be calculated and A and B are constant matrixes.
If we do not consider the item L i 1 ,i 2 ρ n i 1 ,i 2 in Eq. (21) which contains Hamiltonian, d-d interaction and the cooperative dissipation, this equation is simple two-dimension phase diffusion matrix equation. But now Eq. www.nature.com/scientificreports/ (21) contains this item, and it is not a simple phase diffusion matrix equation anymore. That is to say, if we do not consider the first item of Eq. (23), Eq. (23) will actually be Sylvester equation 46 .
The second rearrangement. If  ) T which is a 16 × 1 basis that comes from rearranging , then we will get the separated L i 2 ;i 1 ,i 2 ;i 1 and the separated ρ n i 2 ;i 1 ,1 multiplying: L i 2 ;i 1 ,i 2 ;i 1ρ n i 2 ;i 1 ,1 (column element by column element is also permissible, but the value of the separated L i 2 ;i 1 ,i 2 ;i 1 and the separated ρ n i 2 ;i 1 ,1 will be different from here). Here L i 2 ;i 1 ,i 2 ;i 1 is a 16 × 16 matrix when ρ n i 2 ;i 1 ,1 is rearranged as ρ n i 2 ;i 1 ,1 . Each block, i.e. ρ n i 2 ;i 1 ,1 or ρ n−1 i 2 ;i 1 ,1 , is then 16 × 1 matrix, and the block matrices (ρ n i 2 ;i 1 ,1 ) i 2 ;i 1 ,1 and (ρ n−1 i 2 ;i 1 ,1 ) i 2 ;i 1 ,1 are called vec(X) and vec(C) respectively. Throughout our paper, superscript T means transposition. Under this matrix rearranging, like the way in solving the Sylvester equation, Eq. (25) can be changed into Setting the value A i l = e −i� l = e −i(i l −1) * �� l , A * i l = e i� l = e i(i l −1) * �� l , (l = 1, 2) , Ŵ = γ / , B = −iŴ − g/ , C = −iŴ + g/ , by some calculations on Eq. (26), we can get L i 2 ;i 1 ,i 2  www.nature.com/scientificreports/ Furthermore, following the sequence ( ) i 2 ;i 1 ,1 of the first block rearrangement, the second rearranged block matrix: (L i 2 ;i 1 ,i 2 ;i 1 ) i 2 ;i 1 ,i 2 ;i 1 , is a N 1 N 2 × N 1 N 2 block diagonal matrix, and the diagonal block matrices are L i 2 ;i 1 ,i 2 ;i 1 . So there are N 1 N 2 diagonal blocks in the second rearranged matrix (L i 2 ;i 1 ,i 2 ;i 1 ) i 2 ;i 1 ,i 2 ;i 1 . Thus the first item of Eq. (23), i.e. (L i 2 ;i 1 ,i 2 ;i 1 ρ n i 2 ;i 1 ,i 2 ;i 1 ) i 2 ;i 1 ,i 2 ;i 1 , has been separated into two rearranged matrices (L i 2 ;i 1 ,i 2 ;i 1 ) i 2 ;i 1 ,i 2 ;i 1 and vec(X) multiplexing, and the second item and the last item of Eq. (23) can be rearranged into (I N 2 ⊗A ⊗ I ⊗4 2 ) vec(X) and (B T ⊗I N 1 ⊗ I ⊗4 2 ) vec(X) , respectively. So after this second rearrangement, Eq. (24) can be transformed into where vec(X) or vec(C) comes from rearranging its own N 1 × N 2 square block matrix (i.e., X or C ) column block by column block into N 1 N 2 × 1 column block matrix indexed with subscript i 2 ; i 1 , 1 (i.e., ( ) i 2 ;i 1 ,1 , which is arranged by the sequence value of (i 2 − 1)N 1 + (i 1 − 1) from 0 to N 1 N 2 − 1 just as we have defined), and each block is still 16 × 1 . The detail definitions about the basis of matrices vec(X) , vec(C) and (L i 2 ;i 1 ,i 2 ;i 1 ) i 2 ;i 1 ,i 2 ;i 1 are in the Supplementary Information 1

. Equation (28) can be calculated by the iterations with
This equation can give us the convenience of the iteration calculations. Here the initial value vec(C) = (ρ 0 i 2 ;i 1 ,1 ) i 2 ;i 1 ,1 is N 1 N 2 blocks, and each block is the same. Below we provide each block is ρ 0 i 2 ;i 1 ,1 = (cos 2 θ,0,0, cos θ sin θe −2iφ , 0,0,0,0,0,0,0,0, cos θ sin θ e 2iφ , 0,0, sin 2 θ) T for any 1 ≤ i 2 ≤ N 2 ; 1 ≤ i 1 ≤ N 1 , which is the matrix of the two qubit maximal entangled state |� 0 � = cos θ|00� + sin θ exp(i2φ)|11� . And we provide θ = π/4 and φ = π/5 throughout this paper. Note that this φ is a fixed phase value at the initial time, and it is to be estimated which has no relation with the former phase noise, i.e. 1 and 2 . We can use the program to iterate n m times between vec(C) and vec(X) , i.e., one iteration result vec(X) is used as the initial value vec(C) of the next iteration. In the end we can get the result, a set constituted from a series of vec(X) , namely {(ρ n i 2 ;i 1 ,1 ) i 2 ;i 1 ,1 } (1 n n m ) . Adding the initial state (ρ 0 i 2 ;i 1 ,1 ) i 2 ;i 1 ,1 into this set, we get a new set {(ρ n i 2 ;i 1 ,1 ) i 2 ;i 1 ,1 } for every moment, where n is expand to {0 n n m } . For each n, each ρ n i 2 ;i 1 ,1 , which means with fixed i 2 and i 1 , indicates each evolving result according to its possible initial phase noises 2 and 1 in Hamiltonian of Eq. (9), respectively. Thus for each n, (ρ n i 2 ;i 1 ,1 ) i 2 ;i 1 ,1 has N 1 N 2 independent solutions, or N 1 N 2 block matrices, and each matrix is 16 × 1 . We retransform each of these N 1 N 2 block matrices from 16 × 1 matrix back into 4 × 4 matrix {(ρ n i 2 ;i 1 ,1 ) i 2 ;i 1 ,1 } and take the stochastic average 28,44,45,47 by summing them up and then dividing the summation with N 1 N 2 (these stochastic average is based on that every initial noise possessing any stochastic value in [0, 2π) has equal probability for each qubit). Then we get a series of stochastic averaged density matrixes with discrete form: {ρ n } , (0 n n m ) . ρ n is a 4 × 4 matrix which is the evolved matrix at the moment t. Ultimately, we get the evolved matrix: {ρ n } , (0 n n m ) by our general numeric solution method.
Hence, by using {ρ n } , we can calculate the QFI of two interacting atoms driven by PNLs. By using Eq. and |� n j(k) � are eigenvalue and eigenvector of ρ(φ) n respectively. Note that in Eq. (30a) we need to numerically calculate the derivative of the evolved matrix ρ(φ) n to the estimated phase φ : ∂ρ(φ) n ∂φ . If we change φ of ρ 0 i 2 ;i 1 ,1 into φ + δφ , then using the same way of calculating ρ(φ) n in this paper, we get a evolved matrix ρ(φ + δφ) n , and the In the calculation of this paper, we as well choose δφ = 0.0001 (it www.nature.com/scientificreports/ can be tested mathematically, usually δφ ∼ 0.0001φ , which gives us a better numerical derivative), then we can get the derivative of ρ(φ) n to the estimated phase φ . And subscribe n m + 1 derivatives into Eq. (30a), we can get the QFI of the time discrete form of two atoms driven by the PNLs: {F n Q (φ)} , (0 n n m ).
Discussion of the QFI. Below we use to scale time, i.e., use t to represent t, i.e., we adopt this style in scaling time, which leads to the replacements in above formulas: ω → ω/ , γ → γ / , D → D/ . And for simplicity, we provide the two diffusion coefficients are equal, i.e. D 1 = D 2 = D , and provide N 1 = N 2 = N . And thus we have a 1 = a 2 . Figure 1 indicates QFI and fidelity of two atoms driven by PNLs under the rotating frame, and here we provide D/ = 0.01 , the maximal time t m = 4 when we use the computers to do calculations. And here we do not consider the d-d interaction and the collective decay (i.e., g = γ = 0 ) of the two atoms. And φ = π/5 is the estimated phase (actually QFI has not relation with this value), and the difference interval of the estimated phase is δφ = 0.0001 . We divide [0, 2π) into N uniformly spaced values, i.e., N is the division number of equally dividing [0, 2π) , or the maximal phase noise index in the former derivation of this paper. From the stochastic viewpoint, all values among [0, 2π) have the equal chances to be the initial phase noises, and what is more these phase noises will diffuse with time between these N values {(i − 1) } , i = 1, 2, ..., N . Actually, strictly speaking, 2π needs to be divided into infinite values, which is strict correct, but Limited by the calculation abilities of computers, N can not be too large.
In (a) and (b) of Fig. 1, blue lines are based on the analytical expression of our former work 1 , i.e. the analytical result, and red, green, brown, pink, sky blue lines are the numeric simulation results based on the method of this paper under the same parameters. For the numeric simulation, from Eqs. (15) and (16), we know that when D/ = 0.01 is a provided value, (��) 2 �t needs to be concrete value, which means N 2 /n m will be equal in different simulations. So we can provide N = 34 and n m = 903 (red), N = 28 and n m = 613 (green), N = 25 and n m = 488 (brown), N = 10 and n m = 78 (pink), N = 5 and n m = 20 (sky blue) in Fig. 1a,b. From Fig. 1a, the numerical simulation (red, green, brown, pink or sky blue line) of this paper is close to the analytical result (blue line) of our former work 1 . However, different divided values N reflect the different numeric simulation accuracies on the calculations of the diffusions of two atoms' phase noises. Comparing the red, green, brown, pink, sky blue lines with blue lines in (a) and (b), we find larger N leads to closer results between the numeric simulation result and the analytical result. This is reasonable for us to understand. And this is consistent with the general law of finite element method for numerical solution of physical model, that is, the smaller the spatial step and the time step, the more accurate the simulation results. Because of the phase diffusion rate D/ = 0 , the QFI and the quantum fidelity will be actually a damping oscillation if we consider a long time evolution.
For large interatomic separations k 0 r jk goes to infinity, then g jk = γ jk = 0 for j = k , i.e., there is no coupling between the atoms. And their decay will be independent decay. This case can be solved analytically, as we and other authors have done before 1,32-34 .  In (a,b), blue lines indicate the curve got from the analytical expression, and the red, green, brown, pink, sky blue lines indicate the numerical simulation result got from this paper under the same parameters, respectively. www.nature.com/scientificreports/ When the two atoms are closer, i.e., k 0 r jk is near zero, from the above definitions of g jk and γ jk , we obtain (provide two atomic dipole moments are parallel to each other) 2,49 and the non-diagonal elements of the coefficient matrix of the collective decay can be approximated to be γ , i.e. γ 12 = γ 21 ≈ γ . Below we provide 2 the optical frequency ω ≈ 3 × 10 15 Hz and dipole moment ≈ 1 Debye. In Fig. 2, we provide D/ = 0.01 , N = 20 , n = 400 , g/ = 0 (blue), g/ = 1 (pink), g/ = 2.5 (brown), g/ = 10 (red) and the maximal time is t m = 2 . (a) γ / = 1 4 ; (b) γ / = 1 20 . The estimated phase φ = π/5 , and its difference interval δφ = 0.0001 . Using Eq. (31), we can obtain dipole-dipole coupling of g/γ = 10 corresponds to interatomic distance ≈ 292 nm between two atoms, which is the situation of brown line of Fig. 2a. g/γ = 40 (red line of Fig. 2a), 50 (brown line of Fig. 2b) and 200 (red line of Fig. 2b) correspond to the distance ≈ 184 nm, ≈ 171 nm and ≈ 108 nm, respectively.
From Fig. 2, no mater in (a) or (b), we can find that when other parameters are the same and at the same time. If g ( i.e. the d-d interaction intensity of the two atoms) becomes larger, their QFI becomes larger. This can be used in the open quantum systems to protect the QFI, i.e. at the environment of noises of phase diffusion. This is one result of this paper. Besides, comparing (a) with (b), we can find that smaller decay rate γ leads to larger QFI, when other circumstances are the same. And this is the other result of this paper. Figure 3 is QFI vs. t . We provide g/ = 2.5 , N = 20 , n = 400 , the estimated phase φ = π/5 , its difference interval δφ = 0.0001 , and the maximal time is t m = 2 . Here γ / = 1 20 (red); γ / = 1 4 (brown); γ / = 3/4 (blue). It is easy to obtain, for example, the blue lines of (a) and (b) correspond to interatomic distance ≈ 195 nm between two atoms. From Fig. 3a, it is clear that when the decay rate γ is larger, under the same parameters, the QFI is smaller. So decreasing the decay rate is helpful to increase the QFI of the system. From our former discussion 1 , we find that when 0 ≤ D/ ≤ 2.606 , there exists the information backflow to the system, so the QFI has a recovery process from a small value to a big value with time. Here comparing three curves with different γ in Fig. 3a, the non-Markovian effect can only emerge when the collective decay rate γ is not big, or this effect will be covered up for a big γ . In Fig. 3b, D/ = 5 , according to the former researches, it belongs to Markovian region of the phase noises. Obviously, it does not exist the information backflow and the QFI recovery anymore. And the curves decays all the time, not like the red curve in Fig. 3a which exists the non-Markovian recovery when the collective decay rate is small. www.nature.com/scientificreports/

Conclusion
In this paper, we have studied the dynamical characters of parament estimation of two atoms driven by PNLs.
In the former research 1 we did not consider d-d interaction between the two atoms, and we got the analytical expression of QFI. Here we consider this factor and use the numerical simulation in calculating the QFI under PNLs in this paper. For the circumstance of this PNLs driving, we find that the smaller collective decay rate, the larger QFI of two atoms can be, when other parameters are the same. Besides, when other parameters are the same, the larger d-d interaction of two atoms, the larger QFI of two atoms can be. So we can give a conclusion that collective decay is destructive to QFI of the two atoms, while the d-d interaction is useful in keeping their QFI. This conclusion may be used in the open quantum system to protect QFI under the PNLs environment.
What is more, our scheme gives a general method to solve the model about two atoms with PNLs ultimately, which is suitable for not only two non-interaction atoms but also for two interaction atoms if the Hamiltonian can be written out. And the general numeric solution method of the evolved two-atomic density matrix can be applied to other areas. Knowing the evolved density matrix of the system, many physical problems can be solved, i.e., our general method to calculate the evolved density matrix can be used not only for QFI problems, but also for other physical problems such as fidelity, quantum discord, non-Markovianity and so on.