Quantum coherence [1] arising from the superposition principle is the key factor deviating quantum physics from classical when describing particle systems such as photons [2], which may have a classical wave theory, and matters in classical physics such as atoms [3] or ions [4]. Coherence leads to interference already known in classical optics, and also a purely quantum phenomenon in bipartite and multipartite systems, quantum entanglement [5]. There have been several measures to quantify coherence, and they should fulfill the publicly accepted four conditions of nonnegativity, monotonicity, strong monotonicity, and convexity [1]. One example of coherence measures is the $l_1$ norm of coherence $C_{l_1}$ introduced in Ref. [6]. Note that coherence is dependent on the reference basis we choose. A pure state has zero coherence if it is one of the basis states, but it may have coherence in other bases.

Coherence and entanglement have many similar aspects [7]. For instance, both can be treated as physical resources [1,6,8], and we often need to preserve coherence against decoherence noises as well as entanglement in theories and experiments. The amount of entanglement has been quantified in many ways, such as concurrence [5,9]. Entanglement dynamics in open systems has been studied. Konrad et al. provided a simple relation describing how the entanglement of two-qubit systems evolves as a whole system, under the action of an arbitrary noise channel on one of the components, and the dynamics of the entanglement becomes a very simple form fully captured by an entanglement factorization law [10]. Furthermore, using linear optical setups, Farías et al. experimentally verified the entanglement factorization law under two quite different types of entanglement dynamics [11], and Xu et al. characterized the bipartite entanglement under one-sided open system dynamics [12–14], which includes pure and mixed initial two-photon states under one-sided phase damping and amplitude decay channel [15]. There have been more general results of the entanglement factorization law describing how bipartite high-dimensional entanglement evolves [16], and how multipartite entanglement [17] of a composite quantum system evolves when one of the subsystems undergoes a physical process [18].

Compared with entanglement, the theoretical and experimental results on coherence dynamics in open systems are limited. We need some general law to determine its evolution equation to help us design the effective coherence preservation schemes. Theoretically, Hu et al. introduced a framework of the evolution equation of coherence, and proved a simple factorization relation for the $l_1$ norm of coherence based on this framework, identifying the sets of quantum channels for which this factorization relation holds [19]. The universality of this relation indicates that it applies to many other related coherence and quantum correlation measures [20–22]. Experimentally, the coherence distillation involving strictly incoherent operations (SIOs) [23] and the achievement of the optimal state-conversion probabilities with stochastic IOs [24] have been demonstrated. To our knowledge, experimental results of a general evolution law for coherence have never been reported. Here we provide a factorization law of a definition of coherence for quantum systems under a certain type of quantum operation, and experimentally verify it for qubits and qutrits using a linear optics setup.

An arbitrary pure state in $d$-dimensional Hilbert space $\mathcal{H}$ can be expressed as $|\psi ⟩={\sum }_{i=1}^d{a}_i|i⟩$. When all $a_i=1/\sqrt{d}$, the state is known as the maximally coherent state (MCS) [6] $|{\psi }^{+}⟩={\sum }_{i=1}^d|i⟩/\sqrt{d}$, which should be one of the most coherent states in $\mathcal{H}$. The quantum operation is described by a set of Kraus operators $\{K_n\}$ satisfying ${\sum }_n{K}_n^{†}{K}_n=I$. After the operation, state $\rho$ becomes ${\sum }_n{K}_n\rho K_n^{†}$. If each Kraus operator is diagonal in the reference basis $K_n=\mathrm{diag}(K_{n,1},K_{n,2},\mathrm{...},K_{n,d})$, the operation is known as a genuinely IO (GIO), which has the property that an incoherent state is untouched after acting on it, i.e., for any incoherent state $\delta$, a GIO $\mathrm{\Phi }$ can preserve it $\mathrm{\Phi }[\delta ]=\delta$ [1,25]. Moreover, all Kraus operators of GIO are diagonal in the reference basis [25]. Therefore, the relations between GIO and IO, SIO, maximally IO (MIO) are $\mathrm{GIO}\subset \mathrm{SIO}\subset \mathrm{IO}\subset \mathrm{MIO}$ [25]. In Refs. [26–29], the $G$-concurrence, which is the geometric mean of the Schmidt coefficients for pure states, has been proposed for fully entangled states (i.e., with $d$ nonzero Schmidt coefficients) in $d\times d$ systems. It can describe the entanglement factorization law in bipartite high-dimensional systems [16]. Similarly, one can quantify the full coherence of a $d$-dimensional quantum state $\rho$ by defining the $G$-coherence as $$G(\rho )=d\prod _{i\ne j}{|{\rho }_{ij}|}^{\frac{1}{d(d-1)}},$$(1)which is $d$ times the geometric mean of all $|{\rho }_{ij}|$ with $i\ne j$. When $d=2$ (qubit case), it becomes $2|{\rho }_{12}|$, the same as the $l_1$ norm of coherence $C_{l_1}(\rho )={\sum }_{i\ne j}|{\rho }_{ij}|$ [6]. This definition means the full coherence is zero as long as there exists a zero off-diagonal element in the density matrix, and GIO will not keep the state unchanged if it has nonzero off-diagonal ones. Nevertheless, it will not bring coherence out of nothing. It is worth noticing that neither $G$-concurrence nor $G$-coherence is designed for every entangled state or every coherent state, respectively. Considering a $d\times d$ Hilbert space, only entangled states with $d$ nonzero Schmidt coefficients contain positive $G$-concurrence; other entangled states with less than $d$ nonzero Schmidt coefficients do not contain any $G$-concurrence (such as any entangled two-qubit states in $3\times 3$ Hilbert space). $G$-concurrence is designed for fully entangled states in $d\times d$ systems. Similarly, for a $d$-dimensional system, $G$-coherence is designed for fully coherent states in a $d$-dimensional system.

This is the coherence factorization law: in $d$-dimensional Hilbert space, if a GIO $\mathrm{\Phi }$ acts on an arbitrary state $\rho$, pure or mixed, and produces the final state $\mathrm{\Phi }(\rho )$, then the coherence measure $G$ of the final state is the product of the coherence of the initial state $\rho$ and the MCS after the operation $\mathrm{\Phi }(|{\psi }^{+}⟩⟨{\psi }^{+}|)$, that is, $$G[\mathrm{\Phi }(\rho )]=G(\rho )G[\mathrm{\Phi }(|{\psi }^{+}⟩⟨{\psi }^{+}|)].$$(2)

We can see that the coherence of the final state is defined by two factors, one of which is solely determined by the initial state, and the other is related only to the operation. The operation term is reflected by the coherence of a state, which is similar to Choi–Jamiołkowski isomorphism, also known as channel-state duality [30]. An illustration of the principle is shown in Fig. 1. The proof of the law uses the stronger theorem $${[\mathrm{\Phi }(\rho )]}_{ij}=d{\rho }_{ij}{[\mathrm{\Phi }(|{\psi }^{+}⟩⟨{\psi }^{+}|)]}_{ij}={\rho }_{ij}{[\mathrm{\Phi }(J_d)]}_{ij}$$(3)[or in the Hadamard product form, $\mathrm{\Phi }(\rho )=\rho \circ \mathrm{\Phi }(J_d)$] for GIO for the factorization of each complex-valued matrix element, where $J_d=d|{\psi }^{+}⟩⟨{\psi }^{+}|$ is the $d\times d$ all-ones matrix. The proof of this theorem uses the property of the matrix product that the multiplication of the diagonal matrix $K_n$ on the left means the $i$th row is multiplied by $K_{n,i}$, and that of $K_n^{†}$ on the right means the $j$th column is multiplied by $K_{n,j}^{\ast }$. So the elements of $\mathrm{\Phi }(J_d)$ are $${[\mathrm{\Phi }(J_d)]}_{ij}=\sum _n{(K_n{J}_d{K}_n^{†})}_{ij}=\sum _n{K}_{n,i}{K}_{n,j}^{\ast }.$$(4)

Similarly, $${[\mathrm{\Phi }(\rho )]}_{ij}=\sum _n{(K_n\rho K_n^{†})}_{ij}={\rho }_{ij}\sum _n{K}_{n,i}{K}_{n,j}^{\ast },$$(5)which proves Eq. (3). Multiplying all the off-diagonal elements together and taking the modulus and the $d(d-1)$th root yields Eq. (2).

The coherence of qubit $G(\rho )=2|{\rho }_{12}|$ is irrelevant to the diagonal elements of the density matrix, so quantum state tomography [31] at the $\{|D⟩,|A⟩\}$ and $\{|L⟩,|R⟩\}$ base is enough to yield this value: $$G=\sqrt{{⟨{\sigma }_x⟩}^2+{⟨{\sigma }_y⟩}^2},$$(6)where ${\sigma }_x$ and ${\sigma }_y$ are Pauli operators. In our experiments, we use the path degree of freedom (DOF) of photons as the system, and the polarization is the auxiliary DOF. The initial state is prepared as $|\psi ⟩=\sin 2{\theta }_1|1⟩+\cos 2{\theta }_1|2⟩$, and the Kraus operators of the operation we use are $K_1=\mathrm{diag}(\sin 2{\theta }_2,\cos 2{\theta }_2)$, $K_2=\mathrm{diag}(\cos 2{\theta }_2,\mathrm{i}\sin 2{\theta }_2)$, which are suitable for implementing in an optical setup. When ${\theta }_1=22.5°$, $|\psi ⟩$ is the MCS.

The experimental setup of qubit is shown in Fig. 2(a). An attenuated pulsed 808 nm light from a laser source is collimated into a single-mode fiber and prepared at the horizontal polarization $|H⟩$ by a half-wave plate (HWP) and a polarizing beam splitter (PBS) after emitting from the collimator. Then an HWP with its fast axis at angle ${\theta }_1$ and a beam displacer (BD) produce the initial state $|\psi ⟩$, while the two paths have different polarizations. The HWP at angle ${\theta }_2$ acting on both paths changes the amplitudes at $H$ and $V$ polarizations. A quarter-wave plate (QWP) at 0° is inserted at path 2 to realize the phase factor $i$ in $K_2$. To compensate for the optical path difference, another QWP is inserted at path 1 right after the first BD. Then a thick lithium niobate (${\text{LiNbO}}_3$) crystal introduces an optical path difference between the $H$ and $V$ components larger than the coherence length of the light source from the birefringence effect, and thus destroys the coherence between the two components [32]. To eliminate the polarization difference, an HWP at 22.5° converts $H$ and $V$ polarizations to diagonal ($D$) and anti-diagonal ($A$), respectively. Two HWPs and the second BD select the $H$ component and merge the two paths together, converting the path information into polarization while discarding half the photons. Before merging, two thin glass plates are inserted at the two paths to compensate for the residual phase difference. Then the photons pass a tomography device consisting of a QWP, an HWP, and a PBS, before being coupled into another single-mode fiber, sent to a single-photon detector (SPD), and counted by a computer. Dark counts are subtracted to increase the signal-to-noise ratio.

We measured $G[\mathrm{\Phi }(|\psi ⟩⟨\psi |)]$ with different ${\theta }_1$ and ${\theta }_2$ values using Eq. (6). Then we removed the optical elements for the operation process from the 90° HWP to ${\text{LiNbO}}_3$ and directly performed tomography on the initial states with different ${\theta }_1$ values. The relation between $G$ of the initial state and the theoretical curve $\sin 4{\theta }_1$ is plotted in Fig. 3(a). In Figs. 3(b)–3(e), we plot $G$ of the final states in red dots, the product on the right-hand side of Eq. (2) in blue dots, and the theoretical curves $|\sin 4{\theta }_1\sin 4{\theta }_2|/\sqrt{2}$ with different ${\theta }_1$ values. The product values are close to the directly measured $G$ values, verifying the factorization law.

We use three path DOFs to study the qutrit scenario, where the initial state is $|\psi ⟩=\sqrt{1/3}|1⟩+\sqrt{2/3}\cos 2{\theta }_2|2⟩+\sqrt{2/3}\sin 2{\theta }_2|3⟩$, and we design a phase damping operation between two of the paths and the other one. The corresponding Kraus operators are $K_1=\mathrm{diag}(\cos 2{\theta }_3,\cos 2{\theta }_3,1)$ and $K_2=\mathrm{diag}(\sin 2{\theta }_3,\sin 2{\theta }_3,0)$. Designing a true phase damping operation on all the paths is more complex, as the polarization DOF as an auxiliary is two-dimensional. When ${\theta }_3=0°$, the coherence is preserved, while it is completely destroyed when ${\theta }_3=45°$. The experimental setup is similar to the qubit experiment except for an additional path, a different wave plate setup to realize the operation, and a different tomography method. The essential parts are shown in Fig. 2(b). ${\theta }_2$ and ${\theta }_3$ are the angles of the wave plates controlling the initial state and the operation, respectively. In the tomography process, wave plates, two BDs, and a PBS are needed to project the quantum state to 15 eigenstates of eight Gell-Mann matrices [31,33,34]. The density matrix is calculated from the averages of these Hermitian operators. See Appendix A for more details. Then we can obtain the coherence via $G(\rho )=3\sqrt[3]{|{\rho }_{12}{\rho }_{13}{\rho }_{23}|}$. However, if some of the off-diagonal elements should be zero while others are nonzero, a small value from the experimental errors will cause the measured $G$ to be significantly larger than zero. For example, in an experiment, if the moduli of the three measured off-diagonal elements are 0.01, 0.2, and 0.2, we have $G\approx 0.22$. So, when the theoretical $G=0$, the experimental value will be inaccurate unless most of the off-diagonal elements are zero.

We choose ${\theta }_2=22.5°$ (MCS) and 7.5° (another initial state) and take different ${\theta }_3$ values to perform state tomography. To show the impact of the phase damping operation on the moduli of density matrix elements, we present $|{\rho }_{\mathit{ij}}|$ values from our experiment in Fig. 4(a) when the initial state is MCS and ${\theta }_3=0°$, 15°, 30°, and 45°. Off-diagonal elements ${\rho }_{13}$ and ${\rho }_{23}$ decay as ${\theta }_3$ increases, while ${\rho }_{12}$ and the diagonal elements are roughly unchanged. The coherence values of the initial states are taken when ${\theta }_3=0°$. The measured $G$ values of final states, the products (not applicable for MCS), and the theoretical curves $G(\rho ){|\sin 4{\theta }_2\cos 2{\theta }_3|}^{2/3}$ are plotted in Fig. 4(b). Coherence decays with the increase in ${\theta }_3$. When ${\theta }_3=45°$, coherence should become zero, but the experimental value is inaccurate as we have stated before, and the error calculated from the $3\sigma$ deviation of Poisson distribution at the angle is larger than others.

The law of Eq. (2) still holds for mixed initial states, which are statistical mixtures of different pure states. Weighted averages of count data can be used to simulate the mixed state scenario. We let the initial state be $(1-p)|{\psi }^{+}⟩⟨{\psi }^{+}|+p|{\psi }^{\prime }⟩⟨{\psi }^{\prime }|$ ($0\le p\le 1$), where $|{\psi }^{\prime }⟩$ is the state when ${\theta }_2=7.5°$, and use the weighted averages of count data from $|{\psi }^{+}⟩$ and $|{\psi }^{\prime }⟩$ as the new count data to simulate the mixed state scenario. Under six types of operations where ${\theta }_3=0°$, 7.5°, 15°, 22.5°, 30°, and 37.5°, corresponding to the six solid dots in Fig. 4(b), the calculated $G$ values (light red curves) and the product values (light blue curves) are shown in Fig. 5 when $p$ takes values from zero to one. The two curves are close to each other, verifying the factorization law with mixed input states.

For qudits with higher dimensions, we can still use the path DOF, but the optical setup is more vulnerable to errors from the misalignment of optical devices, making the measured values deviate from theoretical ones. Nevertheless, we performed a four-dimensional experiment, projected the states after operation for qudit state tomography [31], and found the equation still agrees well for a certain operation and initial state: $G[\mathrm{\Phi }(|\psi ⟩⟨\psi |)]\approx 0.4349$, $G(|\psi ⟩⟨\psi |)G[\mathrm{\Phi }(|{\psi }^{+}⟩⟨{\psi }^{+}|)]\approx 0.7406\times 0.5849\approx 0.4332$.

There are other types of operations that satisfy the law. One example is that all the Kraus operators $K_n$ are multiplied by the same permutation matrix on the left, changing the order of the reference bases while keeping the coherence value. For example, a qutrit operation described by $\{K_n^{\prime }\}$ is related to the Kraus operators $\{K_n\}$ of a GIO by $$K_n^{\prime }=\left(\begin{array}{ccc}0& K_{n,2}& 0\\ 0& 0& K_{n,3}\\ K_{n,1}& 0& 0\end{array}\right)=\left(\begin{array}{ccc}0& 1& 0\\ 0& 0& 1\\ 1& 0& 0\end{array}\right)K_n.$$(7)

But it does not hold when the permutation matrices for each $K_n$ are different. Also, some quantum operations cannot be described in the form above, but they are special to satisfy it as well. One example is the qubit amplitude decay channel $$K_1=\left(\begin{array}{cc}1& 0\\ 0& \sqrt{1-ϵ}\end{array}\right),K_2=\left(\begin{array}{cc}0& \sqrt{ϵ}\\ 0& 0\end{array}\right),$$(8)where $0<ϵ<1$, which describes the decay from the excited state to the ground state, whose correctness can be verified through calculation.

In summary, we have presented and proved a factorization law for qudits under GIO, in that using the $G$-coherence measure, the coherence of the state after the operation can be factorized into the product of an initial state term and an operation term. To verify the law, we used an optical setup to test qubit and qutrit cases using a given set of initial states and operations. Our work provides an indirect method to measure the final coherence of quantum states after a specific kind of evolution, and would play an important role in the simplification of coherence measurements, as well as the discovery of other laws about quantum coherence. For example, there are other types of coherence metrics, such as the convex-roof norm (see Ref. [7] for more details) whose calculation method is uncertain yet. We can define $\tilde{G}$ as the minimum statistical average $G$ of all possible pure state combinations of $\rho$: $$\tilde{G}(\rho )=\underset{\{p_i,|{\psi }_i⟩\}}{\min }\sum _i{p}_{i}G(|\psi ⟩⟨\psi |)$$(9){for qubits $\tilde{G}(\rho )=G(\rho )$ [1,7], but calculating this value for qudits is hard}. There may be some factorization law for this coherence measure to be explored in the future, and this would deepen our understanding of quantum coherence. The possible applications of the law may be helping us design the effective coherence preservation schemes in quantum computation and quantum algorithms. For instance, many quantum algorithms have used ${|+⟩}^{\otimes n}$ as input states with $|+⟩=(|0⟩+|1⟩)/\sqrt{2}$, which contain maximal coherence. However, if each qubit of the initial $n$-qubit state passes through a GIO channel such as bit and phase flip together, the output state does not contain maximal coherence. To preserve the coherence of output states, one can use the law and design proper schemes.