The Stimulated Raman Adiabatic Passage (STIRAP), introduced in the 1980s, is a method designed for transferring population between two specific levels in a three-level atomic system [1–4]. STIRAP, along with its expanded theoretical framework, has widespread physical applications [5–9] due to its remarkable properties: (1) during the process of state transfer, the system's evolution remains unaffected by spontaneous emission from intermediate quantum states, allowing for a relatively large decay rate of the intermediate quantum state; (2) it exhibits robustness to minor variations in experimental control parameters, such as pulse amplitude and phase. Among the several extensions of STIRAP, what is of special interest here is fractional STIRAP (f-STIRAP) [4, 10], which involves the simultaneous vanishing of the two pulses (Stokes and pump pulses) with a constant ratio of amplitudes. This approach provides an effective method for creating coherent superpositions of the two levels and can be leveraged to generate entanglement—a fundamental concept in quantum physics and a crucial resource in the development of quantum information and communication technologies [11–14]. Here we rely on f-STIRAP to create a coherent superposition of states $| {e}^{\left(1\right)},{g}^{\left(2\right)}\rangle $ and $| {g}^{\left(1\right)},{e}^{\left(2\right)}\rangle $ of a pair of qubits. Generating entanglement based on adiabatic protocols has been already proposed or successfully demonstrated in various setups, including atomic systems [15, 16], circuit QED [17], NV centers [18, 19], and qubits interacting via a common cavity [20].

In this paper, we consider the optimization of the entangling process in the presence of a realistic dissipative environment. Importantly, we include qubits relaxation besides strong damping of the quantum bus, which makes our analysis especially relevant for solid-state implementations. In these systems, due to the limited qubit life time, dissipation is normally significant for all levels of the Λ-energy structure. Although the intermediate level is most affected, the evolution time in STIRAP processes tends to be lengthy, to satisfy the requirement of being quasi-adiabatic, thus the accumulation of dissipative effects on the ‘dark' state can significantly degrade the desired outcome. Balancing these two competing effects, decoherence and non-adiabatic transitions, becomes crucial for STIRAP under dissipation.

In prior research on the optimization of STIRAP-type protocols, the impact of dissipation has been commonly studied using numerical simulations [21–27]. Here, instead, we adopt the more transparent analytical approach developed in [7], where a STIRAP state transfer process has been analyzed in detail. In that study, the leading-orders corrections induced by weak non-adiabaticity and system's dissipation have been included through a perturbative solution of the quantum master equation. The resulting analytic expressions yield an upper bound for the transfer fidelity, which is in good agreement with direct numerical optimization [7]. With the same approach, the time-dependent coupling profile to realize the optimal state transfer have been obtained as well. Inspired by that work, we apply a perturbative analysis to the f-STIRAP protocol, in order to prepare a maximally entangled state (Bell state) of the two qubits. In this case, the same type of competition between adiabatic and dissipative effects occurs in general, thus requiring us to optimize the quasi-adiabatic entanglement generation protocol under realistic dissipative conditions.

The paper is structured as follows: in section 2 we introduce our system, consisting of two qubits interacting with a bosonic mode, which enables the generation of maximally entangled states between the qubits. In section 3 we discuss the characterization of entanglement between the qubits using concurrence. Section 4 focuses on deriving an approximate expression for the concurrence using a perturbative approach. In sections 5, we optimize the entanglement generation scheme, providing a detailed description of the maximum achievable concurrence and the optimal coupling scheme. Finally, section 6 summarizes our findings and offers an outlook for future research.

We consider two qubits interacting with a common bosonic mode resonantly (see figure 1). In the interaction picture and with the rotating-wave approximation, the Hamiltonian in the interaction picture iswhere G_i(t) denotes the tunable coupling strength between qubit i and the bosonic mode. The operator $\hat{a}$ represents the annihilation operator of the bosonic mode, while ${\hat{\sigma }}_{\pm }^{\left(i\right)}=({\hat{\sigma }}_{x}^{\left(i\right)}\pm {\rm{i}}{\hat{\sigma }}_{y}^{\left(i\right)})/2$ is the spin raising/lowering operator of qubit i (${\hat{\sigma }}_{x,y,z}^{\left(i\right)}$ are Pauli matrices). Such Hamiltonian is widely applied across diverse fields, including quantum optics, quantum information processing, and quantum computing, and can be implemented in various physical systems [28, 29]. The/excitation number ${\hat{N}}_{e}={\hat{a}}^{\dagger }\hat{a}+{\sum }_{i}{\sigma }_{+}^{\left(i\right)}{\sigma }_{-}^{\left(i\right)}$ is conserved by ${\hat{H}}_{I}$ and, similarly to [7], we focus on the low-temperature limit, when excitations cannot be crated by the environment. Based on this, the evolution of system will be confined in a 4-dimensional Hilbert space (see panel (b) of figure 1): $| 1\rangle =| {e}^{\left(1\right)},0,{g}^{\left(2\right)}\rangle ,| 2\rangle =| {g}^{\left(1\right)},1,{g}^{\left(2\right)}\rangle ,| 3\rangle =| {g}^{\left(1\right)},0,{e}^{\left(2\right)}\rangle ,| 4\rangle \,=| {g}^{\left(1\right)},0,{g}^{\left(2\right)}\rangle $. Notice that there is no coherent coupling between state ∣4⟩ and the three other states. In the one-excitation subspace, the instantaneous eigenstates are:with eigenvalues {G(t), 0, − G(t)}, respectively, where $G(t)=\sqrt{{G}_{1}^{2}(t)+{G}_{2}^{2}(t)}$. In equation (2) we used the coupling parameter $\theta (t)=\arctan [{G}_{1}(t)/{G}_{2}(t)]$. Notably, state $| \tilde{2}(t)\rangle $ does not contain the component state ∣2⟩, therefore the bosonic mode has zero occupation. We recall the main feature of STIRAP mentioned in the introduction, i.e., that it allows to avoid dissipation of an intermediate level. Here we shall assume that the bosonic mode suffers a significant decay rate γ. However, state $| \tilde{2}(t)\rangle $ has the desirable property that it neither emits nor absorbs photons to/from the bosonic mode, thus it is referred to as a ‘dark state'. If qubit dissipation mechanisms are negligible, we can adiabatically tune the coupling strength G₂(t) from a finite value to 0, while slowly increasing G₁(t) from 0 to a finite value simultaneously. Doing so, the system will adiabatically evolve from state ∣1⟩ to ∣3⟩ along the dark state $| \tilde{2}(t)\rangle $. In this work, we consider an f-STIRAP process where the parameters are tuned to the middle point G₁(t_f) = G₂(t_f) (giving θ(t_f) = π/4). Under this ideal evolution, the system is prepared in a two-qubit maximally entangled state:

We now include explicitly the effect of environment. Especially important is dissipation of the qubits since, even when it is weak, the long-time adiabatic evolution can lead to significant effects that cannot be ignored. We describe the evolution of the open-system density operator by the quantum master equationwhere $\hat{{ \mathcal L }}={\hat{{ \mathcal L }}}_{m}+{\sum }_{i}{\hat{{ \mathcal L }}}_{q}^{\left(i\right)}$ and the Lindblad superoperators are given by $\hat{{{ \mathcal L }}_{m}}=\gamma { \mathcal D }[\hat{a}]$ and ${\hat{{ \mathcal L }}}_{q}^{\left(i\right)}={\kappa }_{i}{ \mathcal D }[{\hat{\sigma }}_{-}^{\left(i\right)}]$, with ${ \mathcal D }[\hat{A}]\,=\hat{A}\hat{\rho }{\hat{A}}^{\dagger }-\{{\hat{A}}^{\dagger }\hat{A},\hat{\rho }\}/2$.Here, γ (κ_i) is the decay rate of the bosonic mode (qubit i) and we will assume γ ≫ κ_i. To extract the non-adiabatic and dissipative corrections, it will be useful to consider the master equation in the adiabatic basis [7]. By writing $\hat{\rho }^{\prime} (t)={\hat{U}}^{\dagger }(t)\hat{\rho }(t)\hat{U}(t)$, where $\hat{U}={\sum }_{k}| \tilde{k}(t)\rangle \langle k| $ is the time-dependent rotating matrix, equation (4) is transformed to:where $\hat{H}^{\prime} (t)=G(t)(| 1\rangle \langle 1| -| 3\rangle \langle 3| )$ and μ = ∣2⟩⟨1∣ + ∣3⟩⟨2∣ − H. c. . Now equation (5) has clear physical meaning. The Hamiltonian $\hat{H}^{\prime} (t)$ is diagonal, thus state ∣2⟩ (corresponding to the dark state $| \tilde{2}(t)\rangle $ in the original frame) is left unaffected by the dominant part of the unitary evolution. The other two terms describe the the deviations from this ideal evolution. Among them, $\tfrac{\dot{\theta }(t)}{\sqrt{2}}[\mu ,\hat{\rho }^{\prime} (t)]$ introduces transitions between different adiabatic states, caused by finite evolution time in reality. This non-adiabatic term becomes dominant when a fast-rising pulse is applied. On the other hand, the dissipative term ${ \mathcal L }\hat{\rho }^{\prime} (t)$ plays a dominant role if a slow-rising pulse is applied. Both non-adiabaticity and dissipation contribute to the reduction of the population in the maximal entanglement state at the final time t_f. Balancing these two competing effects and optimizing entanglement are the primary goals of our work.

Quantum entanglement plays a crucial role in quantum information processing, particularly in the arenas of quantum communication [11, 12], quantum simulation [13], and quantum computing [14]. Until now, various measures have been proposed to quantify entanglement in different quantum systems [30–34]. Particularly, for the two-qubit quantum system, concurrence [33, 34] has proven to be an efficient entanglement measure. It provides a quantitative metric offering valuable insights into the level of correlation and entanglement present within the quantum state. In our system, the reduced density matrix of the two qubits takes a particularly simple form:where we used the standard basis $| {e}^{\left(1\right)},{e}^{\left(2\right)}\rangle $, $| {e}^{\left(1\right)},{g}^{\left(2\right)}\rangle $, $| {g}^{\left(1\right)},{e}^{\left(2\right)}\rangle $, $| {g}^{\left(1\right)},{g}^{\left(2\right)}\rangle $. The matrix elements ${\rho }_{{ij}}(t)=\langle i| \hat{\rho }(t)| j\rangle $ are given in terms of the states ∣i⟩ defined before equation (2). The form of ${\hat{\rho }}_{R}$ reflects that, as obvious, $| {e}^{\left(1\right)},{e}^{\left(2\right)}\rangle $ does not enter the two-qubit evolution. Furthermore, when tracing out the photon state, coherence between $| {g}^{\left(1\right)},{g}^{\left(2\right)}\rangle $ and the one-excitation states is lost. From equation (2) the concurrence is computed in a standard way as [33, 34]where λ₁, λ₂, λ₃, λ₄ are the square roots of the eigenvalues (taken in decreasing order) of the non-Hermitian matrixIn the above equation (8), ${\hat{\rho }}_{R}^{* }$ is the complex conjugate of ${\hat{\rho }}_{R}$. Finally, we find the simple result:For the maximally entangled state of equation (3), the concurrence is C = 1. Instead, the system is in a product state if C = 0.

Since in the following we will compute the density matrix from the master equation (5), given in the adiabatic basis, it is more appropriate to express ρ₃₁ in terms of the matrix elements of $\hat{\rho }^{\prime} (t)$, by using $\hat{\rho }(t)=\hat{U}(t)\hat{\rho }^{\prime} (t){\hat{U}}^{\dagger }(t)$. We also note that, with the initial condition $| \tilde{2}\rangle $, the two states $| \tilde{1}\rangle $ and $| \tilde{3}\rangle $ play a rather symmetric role in time evolution. For example, their instantaneous eigenvalues ±G only differ by the sign, and the decay rate of $| \tilde{1}\rangle \to | \tilde{4}\rangle $ is exactly the same of $| \tilde{3}\rangle \to | \tilde{4}\rangle $. As a consequence, we find that the matrix elements of $\hat{\rho }^{\prime} (t)$ satisfy ${\rho }_{12}^{{\prime} }(t)=-{\rho }_{23}^{{\prime} }(t)$ and ${\rho }_{11}^{{\prime} }(t)={\rho }_{33}^{{\prime} }(t)$. This allows us to express the concurrence in the following form:From this expression we see that, unlike the state transfer fidelity (given by ${\rho }_{22}^{{\prime} }({t}_{f})$), the concurrence of the final state depends on multiple matrix elements of $\hat{\rho }^{\prime} ({t}_{f})$. Thus, in comparison to the state transfer process, it appears that the optimization of entanglement is more involved. In general, especially when the intermediate state has a significant dissipation rate, it is expected that the evolution is well approximated by $| \tilde{2}(t)\rangle $, thus the parameter θ(t_f) should be close to π/4. In the following, for the analytical solution, we will set θ(t_f) = π/4. Then, the concurrence in equation (10) can be simplified to $C=\left|{\rho }_{22}^{{\prime} }+\mathrm{Re}({\rho }_{13}^{{\prime} })-{\rho }_{11}^{{\prime} }\right|$ when t = t_f.

To evaluate the expression of the concurrence, see equation (10), we utilize the same approach of [7] where, starting from the master equation in equation (5), non-adiabatic and dissipation effects are taken into account perturbatively. We expand the density matrix as $\rho ^{\prime} (t)\,={\rho }^{{\prime} (0)}(t)+{\rho }^{{\prime} (1)}(t)+{\rho }^{{\prime} (2)}(t)+...$, where the zero-order term satisfiesand the matrix element of the higher-order terms satisfy:where we defined E_ab(t) = E_a(t) − E_b(t), with E_a(t) = G(t)(δ_a,1 − δ_a,3). In the integrand of equation (12), ${\xi }^{\left(k-1\right)}(\tau )=[\mu ,{\rho }^{{\prime} (k-1)}(\tau )]$ describes transitions between different adiabatic eigenstates. Instead, the dissipation term is ${L}_{{ab}}^{\left(k-1\right)}(\tau )={\left({ \mathcal L }{\rho }^{{\prime} (k-1)}(\tau )\right)}_{{ab}}$. Based on this perturbative expansion, it will be necessary to obtain explicit expressions of the density matrix up to the third-order. For simplicity we will assume θ(t_f) = π/4, thus we only need to consider the matrix elements ${\rho }_{11}^{{\prime} }(t)$, ${\rho }_{13}^{{\prime} }(t)$, and ${\rho }_{22}^{{\prime} }(t)$.

At zero-order (k = 0), all non-adiabatic and dissipation effects are disregarded and throughout the entire evolution process the system remains in the dark state $| \tilde{2}\rangle $. Consequently, we have ${\rho }_{22}^{{\prime} (0)}=1$, while all other matrix elements are zero. In particular, ${\rho }_{11}^{{\prime} }(t)={\rho }_{13}^{{\prime} }(t)=0$. Instead, the expressions for the first-order (k = 1) density matrix should be derived from equation (12). In particular, we obtain:Iterating again equation (12), the second-order perturbation (k = 2) can be obtained. Typically, we assume that dissipation of the qubit is relatively small, thus its effect can be neglected in higher-order perturbations if they have already been accounted for in ${\rho }^{{\prime} (1)}$. For the relevant matrix elements we find:Finally, the third-order perturbation (k = 3) can be derived using the same method. The resulting expressions have a structure similar to equation (14) but are significantly more involved and not particularly instructive, so we omit them. A simplified form, useful to calculate the concurrence, will be given in equation (17) below.

We now simplify the expressions of the various perturbative terms (k = 1, 2, 3), which comprise integrals of the form ${\int }_{0}^{t}{\rm{d}}t^{\prime} f(t^{\prime} )$ $\exp [\pm {\rm{i}}{\int }_{0}^{t^{\prime} }{\rm{d}}\tau G(\tau )]$. Within these integrals, $f(t^{\prime} )$ is regarded as a relatively smooth function, while $\exp [\pm {\rm{i}}{\int }_{0}^{t^{\prime} }{\rm{d}}\tau G(\tau )]$ is a rapidly oscillating term. This is due to our assumption of adiabaticity, requiring the dynamics of the system to be much slower that the fast timescale associated with the large coupling strength G(t). Therefore, we can simplify these expressions by performing integrations by parts and disregarding the higher-order terms, suppressed by G(t)⁻¹. At the final time t = t_f we get:while the second-order expressions are:We finally give the simplified form of the third-order expressions:In closing this section, we note that the expressions given above for ${\rho }_{22}^{{\prime} }$ coincide with the ones in [7]. Instead, the results for ${\rho }_{11}^{{\prime} }$ and ${\rho }_{13}^{{\prime} }$, necessary to compute the concurrence, are new.

Based on the perturbative expansion of the density matrix, we consider the optimization of entanglement for the case of a fixed final mixing angle θ(t_f) = π/4. In general, one of the most common coupling schemes in STIRAP is the parallel adiabatic passage (PAP), where ${G}_{1}(t)={G}_{0}\sin \theta (t),{G}_{2}(t)\,={G}_{0}\cos \theta (t)$. For this choice of of couplings, the energy gap is constant (G(t) = G₀), which helps in suppressing non-adiabatic errors by preventing energy-level crossings. Furthermore, to achieve the level of high-fidelity operations considered here, it is reasonable to expect that the qubits are obtained from a controlled fabrication process, leading to very similar parameters. We will assume identical decay rates, κ₁ = κ₂ = κ, which also allows us to obtain simple analytical expressions for the optimal concurrence. Under these assumptions, we find:where ${A}_{1}=(2-\gamma {t}_{f}/2)\cos ({G}_{0}{t}_{f})$ and ${A}_{2}=(-1+\gamma {t}_{f}/2){\sin }^{2}({G}_{0}{t}_{f})$. This expression for the concurrence is very similar to the fidelity of [7]. In particular, the important terms in parenthesis arises from ${\rho }_{22}^{{\prime} }$, thus they also appears in the expression of the fidelity. Instead, the remaining boundary terms include contributions of ${\rho }_{11}^{{\prime} }$ and ${\rho }_{13}^{{\prime} }$, absent in the fidelity. In both cases, the concurrence is only related to the time-dependent $\dot{\theta }(t)$, which we can express through a truncated Fourier series:Above, to satisfy the boundary condition θ(t_f) = π/4, we should fix α₀ = π/4t_f. Instead, the other α_n is optimizing parameters which, by imposing ∂C(ρ)/∂α_n = 0, can be obtained explicitly:As seen, all the odd-order and all the even-order are the same, respectively, and the even-order terms dominate at large N. Inserting these coefficients into equations (19) and (18), we find the optimized concurrence:Finally, we consider the limit of a Fourier expansion with N → ∞ , which accounts for a $\dot{\theta }(t)$ with arbitrary time-dependence, thus corresponds to a full optimization of the concurrence:The two terms reducing the concurrence clearly reflect the competition between qubit dissipation (∝κt_f) and non-adiabaticity (∝γ/t_f). We can further balance these two effects by selecting the optimal evolution time ${t}_{f}^{\mathrm{opt}}=\pi /\left(2\kappa \sqrt{4{G}_{0}^{2}/\gamma \kappa }\right)$, which gives:The upper bound of concurrence is closely related to the cooperativity $4{G}_{0}^{2}/\gamma \kappa $. We note that the above expressions differ from the corresponding formulas of the fidelity only by simple numerical factors 1/4 and 1/2, which multiply the non-adiabatic term of equation (22) and the total loss of concurrency of equation (23), respectively. These differences are due to the boundary condition, θ(t_f) = π/4, instead of θ(t_f) = π/2.

In figure 2(a), we present a comparison between optimized numerical and analytical solutions. The blue, purple, magenta, and green curves correspond to numerical optimizations performed at different orders N = 0, 1, 2, 4 of the Fourier expansion, i.e., including the 2N coefficients {α₁, α₂, …α_2N}. The black dashed line represents the analytical solution for optimized entanglement, see equation (22). The horizontal red line indicates the analytical upper limit for the concurrence given by equation (23). From this figure, we can conclude that perturbation theory is unsuitable when the adiabatic evolution time t_f is relatively short, due to significant non-adiabatic effects. Similarly, when the evolution time t_f is large, dissipative effects become prominent, also rendering perturbation theory inapplicable. However, our analytical and numerical solutions are in good agreement in the intermediate range, where the largest values of the concurrence are attained. The numerically obtained time-dependence of θ(t), with t_f close to the optimal condition (G₀t_f ≈ 8), is depicted in panel figure 2(b).

We can further optimize the entanglement numerically by considering the value θ(t_f) as an additional optimization parameters, i.e., we relax the condition θ(t_f) = π/4 assumed so far. From the point of view of the Fourier expansion (19), this is equivalent to include α₀ in the parameters to optimize. As shown in figure 3(a), however, the gain of concurrence is modest. One can see that the difference between the concurrence with θ(t_f) ≠ π/4 (red dots) and θ(t_f) = π/4 (blue dots) is always relatively small, and gets progressively reduced as the expansion order N increases.

Despite the fact that for given parameters the two series of values (blue vs red) appear to give nearly identical results when N → ∞. We find that, surprisingly, the numerically optimized values of θ(t_f) significantly differ from π/4. As seen in figure 3(b), however, the pulse shape does not really deviate strongly from the θ(t_f) = π/4 curve. Considering the explicit form of θ(t), we see that it indeed gradually approaches π/4 and a strong departure only occurs in a short time interval before t_f. For the red curve, differently from the optimized dependence with a fixed θ(t_f) = π/4 (blue curve), there is a sudden decrease in amplitude at the end of the pulse, which is responsible for the slight enhancement of cooperativity. Such an effect, induced by a strong deviation from adiabaticity, cannot be captured by our analytical treatment. However, in the limit of a small decay rate of the qubits, it is expected that the adiabatic evolution becomes nearly optimal, and the dominant loss of cooperativity occurs in the bulk of the time interval 0 < t < t_f, rather than close to t_f. Therefore, the main contributions to the loss of cooperativity should be accurately captured by our treatment. The numerical results depicted in figure 3 confirm these arguments, showing that the appropriate parameter regime is not necessary to account for such non-adiabatic effects in the analytical treatment and experimental realization, as their impact is marginal.

We have explored the generation of maximally entangled states of two qubits interacting with a common quantum bus, based on an adaptation of f-STIRAP to the present setup. We analytically derived an approximate expression for concurrence by employing perturbation theory. Additionally, we optimized the entanglement by balancing non-adiabatic and dissipative effects during the generation of entangled states, and provided the corresponding optimal coupling strategy. This work can serve as theoretical guidance for efficiently realizing two-qubit maximal entangled states based on f-STIRAP, especially for solid-state systems in which qubit dissipation might have a significant effect during the course of the time evolution. Additionally, it can also be principally extended to multipartite entanglement generation. In future work, we aim to extend our treatment using the Lagrangian formalism method of [35], which can be applied to the optimization of quasi-adiabatic evolution under a more general dissipative environment.