Fast quantum state transfer and entanglement preparation in strongly coupled bosonic systems

Continuous U(1) gauge symmetry, which guarantees the conservation of total excitations in linear bosonic systems, will be broken when it comes to the strong-coupling regime where the rotation wave approximation (RWA) fails. Here we develop analytic solutions for multi-mode bosonic systems with XX-type couplings beyond RWA, and propose a novel scheme to implement high-fidelity quantum state transfer (QST) and entanglement preparation (EP) with high speed. The scheme can be realized with designated coupling strength and pulse duration with which the excitation number keeps unchanged regardless of the breakdown of the global U(1) symmetry. In QST tasks, we consider several typical quantum states and demonstrate that this method is robust against thermal noise and imperfections of experimental sequence. In EP tasks, the scheme is successfully implemented for the preparation of Bell states and W-type states, within a shortest preparation time.


Introduction
In the past several decades, the linear transformation of boson modes has been an intriguing topic in quantum optics and quantum information sciences.For example, it plays a crucial role in the problem of Bose samplings [1][2][3][4][5], which is of particular interest in the study of near-term platform for photonic quantum computing [6].The linear bosonic transformation can be effectively realized by optical beam splitters and wave plates in optical systems under the protocol raised by Knill, Laflamme and Milburn [7][8][9], or by controllable pulse manipulation in superconductor-waveguide systems [10][11][12][13][14][15].Experimental progresses in this direction have facilitated the development of various quantum information tasks, including quantum state transfer (QST), entanglement preparation (EP) and entanglement distribution [12][13][14][15][16][17].
As an important step of quantum information processing, QST aims to transfer an arbitrary quantum state from the sender side to the receiver side with high fidelity and fast speed.Recently, to embrace the noisy intermediate-scale quantum era for quantum internet frameworks [15,[18][19][20][21], much effort has been made to implement QST tasks in various physical systems, including atom-cavity systems [10,[21][22][23], superconducting circuits [10][11][12][13][14]24], photonic systems [25][26][27], mechanical oscillators [28], optomechanical cavities [16,[29][30][31][32][33] and so on [17,[34][35][36][37].Meanwhile, the generation of entangled states, as the first and fundamental step to realize quantum algorithms and manifest the so-called quantum supremacy, has been widely studied [12,14,[38][39][40][41].By manipulating the bosonic transformation matrix, entanglement resource can be generated and distributed among different ports of user.In the weak-coupling regime where the rotation wave approximation (RWA) can be safely adopted, QST and EP can be accomplished with high fidelity but slow speed [10,[12][13][14].On the other hand, as one tries to increase the processing speed and pushes to the strong-coupling regime, the breakdown of the global U(1) symmetry will lead to the failure of RWA, resulting in the deviation of the desired bosonic transformation.In such a case, one has to go beyond RWA and develop new schemes to suppress errors brought by the broken U(1) symmetry.
In this paper, we develop a scheme to realize linear bosonic transformation for QST and EP in the strong-coupling regime.By formulating an analytic solution of the widelyused multi-mode model, we find that the total excitation number will be preserved for accurately controlled coupling strength and at certain discrete points of time, in despite of the broken U(1) symmetry.For QST tasks involving two terminal modes and one intermediate channel mode, we obtain an analytical tradeoff relation between the transfer speed and the fidelity.Taking several typical states as examples, we further demonstrate that our method outperforms the traditional one derived from RWA with lower infidelity and absolute robustness against thermal noise of the intermediate mode.The fidelity can be further improved by applying a simple local operation to compensate the phase rotation induced by the strong-coupling.And we also show the potential to extend our protocol to multi-mode QST tasks through an example of transferring an arbitrary W-type state.For EP tasks, we show that typical entangled states such as Bell states and multi-mode W-type states [42] can be successfully generated using the proposed scheme, where the fastest preparation time with our scheme can also be derived.The degree of entanglement can be modified by changing the coupling strengths of different modes.

Model
We start from a widely adopted model where n oscillator-encoded qubits are coupled to an intermediate channel mode with an XX-type coupling [28,29,43,44].This model is widely applied to implement QST protocol in various setups such as opto-mechanical systems and macroscopic harmonic oscillators [45,46].The Hamiltonian can be written in the interaction frame as (we set = k B = 1 in this paper) The boson modes a j represents the jth oscillator-encoded qubits around the common central mode c.It's also assumed that all the modes are resonant with frequency ω.The coupling strengths between the boson modes and the central mode can be effectively controlled by rectangle pulses focused on different mode, and the relative amplitude of the jth pulse is denoted by k j .
The multi-mode bosonic model can be realized in circuit quantum electrodynamical systems, such as superconductor chips [12,13], where spatially separated superconducting qubits are connected by inductor-capacitor (LC) circuit [47].The LC circuit can serve as a single-mode quantum data bus with high quality factor for transferring quantum information between superconductor qubits.In addition, every superconductor qubit individually couples to a single-mode readout cavity.Thus, by connecting the n readout cavities with a common LC circuit one can realize the model under consideration.
Distinguished from other approaches based on RWA in the long-time limit t ≫ 1/ω, the counterrotation terms a j ce −2iωt and a † j c † e 2iωt are retained in our model, thus breaking the global U(1) symmetry.Such generalization makes it a more realistic model in the strong-coupling regime, or equivalently, short-time limit t ∼ 1/ω.By solving the Heisenberg equation of the system, we can obtain an analytical expression of the linear transformation for all the boson modes as a function of evolution time t (see Appendix A).
In the following sections, we will first focus on a simple case of n = 2 to demonstrate our designed fast linear bosonic transformation.In such three-mode (two boson modes plus one channel mode) system, the analytical solution can be obtained by defining a where ] is determined by the initial condition.The solution of mode c can be directly obtained according to the equivalent position between c(t) and the united mode a(t).

QST of a single boson mode
In this section, we discuss the application of the analytic solution to tasks of transferring a single boson mode from the sender a 1 to the receiver a 2 through the intermediate channel mode c, i.e., taking the simplest example of n = 2 in Eq. (1).For simplicity, we consider equal couplings k 1 = k 2 = 1, while the qualitative conclusions can be generalized to other cases.In the weak-coupling regime, the implementation of QST has been studied using conventional method based on RWA [10,13,17,23,27,32,35].By neglecting the counterrotation terms a j ce −2iωt and a † j c † e 2iωt in the long-time limit t ≫ 1/ω, the initial Hamiltonian is reduced to Hint = 2 j=1 g(a j c † + ca † j ), with global U(1) symmetry and conserved total excitation number.Then, the Heisenberg equations can be obtained as Here the coupling strength g and the duration τ of the rectangle pulse should satisfy the relation g ′ τ = π, where g ′ ≡ √ 2g is the effective coupling strength determined by the pulse amplitude [13,17,23,27,32,35].Under such condition, the exchange of modes a 1 and a 2 can be realized, a 1 (τ ) = −a 2 (0), a 2 (τ ) = −a 1 (0), while the intermediate mode c remaining unchanged, c(τ ) = −c(0).However, this perfect result of QST is only an artifact rooted from the assumed global U(1) symmetry which is only approximately preserved for weak coupling.

Optimized bosonic transformation for fast QST
When aiming towards high-speed QST, a strong coupling is required and counterrotation terms become non-negligible, which break the global U(1) gauge invariance and the conservation of total excitation number, and lead to sizable infidelity of the final state.Thus, one needs to go beyond RWA and work with the original Hamiltonian.
Defining the vector of the operators as Ψ = (a 1 , c, a 2 ) T and T , the dynamics of all the three modes right after a duration t can be expressed by the block matrix form in the Heisenberg picture, Typically, the creation and annihilation vectors couple with each other by the offdiagonal transformation matrices U B (t) and U * B (t), and cause the variation of total excitation number.We focus on the first line of the 6-dimensional linear expression of Eq. ( 5), and write the final state of mode a 1 right after the pulse duration τ .
where the coefficients K's can be expressed as functions of τ and g ′ (see Appendix B for details).A successful QST with high fidelity corresponds to the case with |K 21 | ∼ 1 and all other coefficients |K| ∼ 0. At the same time, in despite of the breakdown of U(1) gauge invariance, we find that the excitation number can be preserved in specific conditions with K 12 = K 22 = K c2 = 0.The optimized scheme of QST can be obtained by solving these equations, leading to Here, m = 2, 3, . . . is an arbitrary positive integer.One can easily find that ζ and θ increase monotonously with parameter m.In (b) and (c), the integer m is taken from 5 to 17, corresponding to an evolution time τ from 2× 2π ω to 8× 2π ω .The solid dots and dashed lines stand for the optimized method and RWA method, respectively.The analytical prediction of the error Eq. ( 11) agrees well with the solid dots, being indistinguishable in this plot (see Appendix D).
The condition Eq. ( 7) provides an implicit constraint to optimize the QST procedure.First, one needs to choose a suitable parameter m.The amplitude and duration of the rectangle pulse are then determined by ζ and θ, respectively, both acquiring discrete values.As shown in Fig. 1(a), the optimized pulse parameters agree well with those obtained within RWA in the weak-coupling limit of small g ′ , but deviate in the strong-coupling regime.Under the optimized pulse condition, the final mode can be simplified as a 1 (τ ) = K 11 a 1 (0) + K 21 a 2 (0) (see Appendix B), where This result suggests that both the amplitude error and phase error of QST are dependent on a single parameter, the phase shift θ r , which is analytically determined by the pulse parameters and can be compensated by a local rotation as we discuss latter.With this pulse ansatz, the annihilation and creation operators are decoupled temporarily, i.e.U B (τ ) = 0, indicating that the total excitation number is conserved at this point although the global U(1) symmetry is broken.A detailed discussion is provided in Appendix C. The linear transformation of the bosonic operators takes the form with transform matrix making it possible for high-fidelity QST.
To demonstrate the performance of the optimized scheme, we consider a QST task from the initial state |Ψ(0) = |ψ, 0, 0 to the goal state |Ψ goal = |0, 0, ψ .The three parts in Dirac kets represent the state in node a 1 (sender), channel c and node a 2 (receiver), respectively.As an example, we first demonstrate the results of transferring a Fock state |ψ = |n with n = 1, 2, 3 [Fig.1(b)], and transferring a coherent state |ψ = |α with α = 0.6e iπ/2 , e iπ/2 , 1.4e iπ/2 [Fig.1(c)] through an ideal channel at zero temperature.By choosing the parameter m from 5 to 17, the infidelity of the received state 1 − f changing with the pulse duration τ is presented, where the fidelity f ≡ tr(  [48].The optimized pulse given by Eq. ( 7) outperforms the one predicted by RWA for all cases.In the weak-coupling regime g ′ 0.1ω, where the parameter is chosen as m 11 and the pulse duration τ 5 × (2π/ω), the optimized scheme reduces to the pulse ansatz under RWA.In the strong-coupling regime g ′ 0.1ω with m 11, our scheme significantly suppresses infidelity from the RWA result, which presents significant fluctuation owing to the non-negligible effect of the counterrotation terms.The reduction of infidelity over the best performance of RWA can be as large as ∼20%, which is substantial considering the already low enough baseline.The ultra strong-coupling regime with m = 2, 3, 4 is not shown in Fig. 1 since the infidelity will be too large to qualify a successful QST process.
Another advantage of our optimized scheme is the independence on the phase of initial state when transferring coherent states.For an initial state ρ = |α α| with α = |α| e iϕ , the infidelity will not change with the phase ϕ as shown in the inset of Fig. 1(c).As a comparison, for the RWA scheme, the infidelity oscillates with ϕ.This is because in RWA method, the final mode a † 1 (τ ) contains not only the target mode a † 2 (0) but also a portion of mode a 2 (0).The additional mode will result in a state similar to the target coherent state |α but with a loss of two photons in every Fock basis, which is expressed as The coefficient iE is pure imaginary with fixed phase π/2, and the real coefficient ε n in the summation is obtained from the binomial expansion [K 21 a † 2 (0) + K 22 a 2 (0) + . ..] n .By combining the effects of these two modes, a relative phase of δθ = 2ϕ + π/2 will emerge and contribute a periodic modulation of infidelity.As a comparison, in our optimized method with the conditions (7) satisfied, the final mode a † 1 (τ ) only contains the target mode a † 2 (0), such that the infidelity is not affected by the phase ϕ.

Tradeoff between speed and fidelity
One key merit of having the analytic solution Eq. ( 8) is that the explicit relation between transfer speed and fidelity of QST can be obtained for the optimized scheme of rectangle pulses.We derive a general constraint of the initial pure state ρ to be transferred, the To demonstrate the performance of the optimized scheme, we consider a QST task from the initial state |Ψ(0) = |ψ, 0, 0 to the goal state |Ψ goal = |0, 0, ψ .The three parts in Dirac kets represent the state in node a 1 (sender), the channel c and the node a 2 (receiver), respectively.As an example, we study the cases of transferring the Fock state |ψ = |n with n = 1, 2, 3 [Fig.1(b)], and transferring the coherent state |ψ = |α with α = 0.6e i π 2 , e i π 2 , 1.4e i π 2 [Fig.1(c)] through an ideal channel at zero temperature.By choosing the parameter k from 5 to 17, the infidelity of the received state 1 − f changing with the pulse duration τ is presented, where the fidelity f ≡ tr( ρ i = |ψ ψ| the density matrix of the sending state at node a 1 and ρ f ≡ tr 1,c (ρ 1,c,2 (τ )) the reduced density matrix of the received state at node a 2 [39].Notice that the optimized method outperforms the RWA method for all cases considered here.In the weak-coupling regime g ′ 0.1ω, where the parameter is chosen as k 11 and the pulse duration τ 5 × 2π ω , the optimized method reduces to the best performed pulse ansatz of RWA method.In the strong-coupling regime g ′ 0.1ω with k 11, the infidelity is significantly suppressed comparing to the RWA method.The ultrastrong-coupling regime with k = 2, 3, 4 is not shown in Fig. 1 since the infidelity will be too large to quantify a successful QST process.
Error-dependent QSL.It is of great importance to quantify the tradeoff between the transfer speed and the error emerged from the environment noise and system defection in QST tasks.For the optimized scheme of rectangle pulses, we obtained a general constraint of the initial pure state ρ to be transferred, the evolution time τ , and the infidelity 1 − f , in the form of where the parameter k is taken as a continuou temporarily.Thus, for a task of sending state upper bound of tolerable error the threshold of θ can be obtained by θ th = θ(⌊ where ⌊k th ⌋ is the nearest integer less than or k th .The QSL, i.e., the fastest possible time plish the transfer, is τ th = θ th /ω.We demon detailed comparison between the infidelity in mized method and that predicted by the error-d QSL as Fig. 1 in SM [34], which coincide well other.
Robustness analysis.To demonstrate the fea the proposed scheme under realistic experimen tions, next we investigate the effects of therma the channel, imperfection of pulse shape, and th phase of the initial state.
For the thermal noise, suppose the initia the channel is a thermal state ρ c = 1 Z e −β than a vacuum state at zero temperature, w 1/T and Z = tr(e −βH c ) is the partition fu Maxwell-Boltzmann distribution.Here we set mann constant k B = 1 for simplification.In Fig compare the optimized method (red squares an angles) with the RWA method (red and blue das at finite temperatures.For the RWA scheme delity increases almost linearly with temperatu ing a worse performance when transferring a F than a coherent state.As a comparison, our evolution time τ , and the infidelity 1 − f , in the form of According to Eq. ( 8), we can define the magnitude of coefficient where m is taken as a continuous variable temporarily.For a task of sending state ρ with an upper bound of tolerable error E tol ≥ 1 − f , the QST time τ must satisfy G(m) ≤ E tol / n 1 ρ .By solving the inverse function of G(m), we can get the threshold (lower bound of m) m th ≡ G −1 ( E tol / n 1 ρ ).Thus, the threshold of θ can be obtained by θ th = θ(⌊m th ⌋ + 1), where ⌊m th ⌋ is the nearest integer less than or equal to m th .The fastest possible time to accomplish the transfer, i.e., the quantum speed limit of QST, is τ th = θ th /ω.We demonstrate a detailed comparison between the infidelity in our optimized method and that predicted by the tradeoff relation in Fig. D1 in Appendix D, which shows excellent agreement for all parameters.

Robustness against fluctuations
To demonstrate the feasibility of the proposed scheme under realistic experimental conditions, in the following we investigate the thermal noise effect of the channel, and the imperfection of pulse shape.
For the thermal noise, suppose the initial state of the channel is a thermal state ρ c = e −βHc /Z rather than a vacuum state at zero temperature, where β = 1/T is the inverse temperature, Z = tr(e −βHc ) is the partition function of Maxwell-Boltzmann distribution, and H c is the channel mode Hamiltonian.In Fig. 2(a), we compare the optimized scheme with RWA method for temperature up to T ∼ 3/ω.While the infidelity of the RWA scheme increases almost linearly with temperature, and showing a worse performance when transferring a Fock state than a coherent state, the optimized method is completely immune to thermal noise in both cases.To understand this observation, we consider as an example the task of transferring an initial Fock state |n 1 through a channel with an initial phonon number n c .Note that at the end of the pulse duration τ , the final state of the nodes are completely decoupled from the channel, leading to a 1 (τ ) = K 11 a 1 (0) + K 21 a 2 (0), a 2 (τ ) = K 11 a 2 (0) + K 21 a 1 (0), and c(τ ) = −e 2iθr c(0).The final state is given by Here, U(τ ) is the unitary operator describing the evolution of the system with time τ .This result shows that after the pulse is applied, the channel mode evolves back to its initial state, being independent on the status of nodes a 1,2 (0).Meanwhile, the influence of the channel on the nodes is also erased at this exact time after a partial trace.This conclusion can be easily generalized to a thermal channel, which is a classical superposition of Fock states with different n c , and also to the case of transferring a coherent state as it can be expanded into Fock basis.
For fluctuation effect induced by the pulse, we stress that in the more interesting strong-coupling regime with large pulse amplitude and short duration, a small deviation of coupling intensity ∆g should be less significant in comparison to an error of duration ∆τ .In Fig. 2(b), we impose different amount of ∆τ , and plot the most prominent infidelity 1 − f for a pulse duration within [τ − ∆τ, τ + ∆τ ].Here, the fidelities for different |∆τ | are displayed by symbols, which are connected by solid lines as a guide for the eyes.As expected, the fluctuation of pulse time is more influential for stronger coupling and shorter transfer time.However, if the fluctuation is relatively small with ∆τ ≤ 0.05×2π/ω, our optimized method still performs well and the increase of infidelity is restricted within one percent even for the strongest coupling considered.

Correction of the phase error
Based on the analytic bosonic transformation, we can obtain the phase rotation of the final mode as in Eq. ( 8).This effect is more severe in the strong-coupling regime, which induces a sizable error to the final state.With that knowledge, one can further enhance the fidelity of QST by applying a local rotation to node a 2 .It can be realized by implementing a local pulse on node a 2 after the coupling pulse (7), which is expressed as H r = g r (t)a † 2 a 2 with the pulse area g r (t)dt = θ r .Then the corrected bosonic vector and the real coefficient ε n in the summation is obtained from the binomial expansion [K 21 a † 2 (0)+K 22 a 2 (0)+. . .] n .By combining the effects of these two modes, a relative phase of δθ = 2ϕ + π/2 will emerge and contribute a periodic modulation of infidelity.However, in our optimized method with the conditions (3) satisfied, the final mode a † 1 (τ ) only contains the target mode a † 2 (0).Hence, the infidelity is not affected by the phase of the initial coherent state ϕ.
Correction of the phase errors.Based on the analytic solution of the equation of motion, we can obtain the phase rotation of the final mode as in Eq. ( 4).This effect is more severe in the strong-coupling regime, and can induce a sizable error to the final state.With that knowledge, one can further enhance the fidelity of QST by applying a local rotation to node a 2 .It can be real-  E1 in Appendix E, such a simple local rotation improves the transferring fidelity remarkably for a coherent state.However, this method does not make any difference for a Fock state since the Wigner function of which is central symmetric.
To elucidate it more explicitly, we consider a task of transferring an even cat state in continuous variable systems |cat(α) + ≡ |α + |−α through an ideal zero temperature channel.As a typical macroscopic quantum superposition, cat state has been prepared in various physical platforms including optical systems [49][50][51], optomechanical systems [52], superconducting systems [53], atomic ensembles [54,55] and magnon-photon systems [56].In Fig. 3(a), we apply our optimized method directly and then observe the final state in node a 2 , which presenting a rotation of an extra phase as expected.In Fig. 3(b), we apply the local rotation to node a 2 right after the coupling pulse (7), and the fidelity is obviously increased.

QST for multi-mode W-type state
Next, we show that the protocol developed for single-mode QST can be extended for transferring multi-mode quantum states.For definiteness, we consider as an example an n s -mode W-type state, which is denoted by |ψ s = ns j=1 C j |1 j ⊗ |0j with |0j ≡ ⊗ i =j |0 i .Without loss of generality, we assume all coefficients C j (j = 1, 2, . . ., n s ) are real, since a complex amplitude can be realized by applying corresponding local phase gate or equivalently by a proper redefinition of the corresponding boson mode a j → a j e iφ .The initial state of the whole system can be expressed as |ψ(0) sr = |ψ s ⊗ |0 r , where the receiver side is assumed to be an n s -mode vacuum state.In the following, we denote the n s boson modes in the sender side as a 1 , . . ., a ns , and the n s boson modes in the receiver side as a ns+1 , . . ., a 2ns to simplify notation.Notice that since the initial state of the channel mode |ψ c will have no influence on the final state, we can trace out this degree of freedom and focus on the 2n s boson modes in sender and receives sides only.Using this notation, the initial and final states can be obtained as where |0 sr = |0 s ⊗ |0 r denotes the vacuum state of both the sender and receiver.By using the complete expression in Appendix A.2, an exact result of the final state can be written down and an optimized transfer scheme can be obtained.On the other hand, to demonstrate the feasibility and advantage of our proposed scheme, in the following discussion we ignore the error caused by the strong coupling temporarily, and write the approximate linear bosonic transformation for the 2n s modes system as Combining with the expression of Eq. ( 14), the final state can be expanded with the Fock basis of initial boson modes, leading to The analytic expressions of the coefficients of the 2n s modes thus allow us to realized QST by properly designing the coupling strengths k j 's.Firstly, we demand that excitations should not occur in the sender side.This requirement gives a total of n s restrictions for the n s coefficients of the sender modes where j 0 = 1, 2, . . ., n s .Then, the relative coupling strengths in the receiver side should be designed in a certain proportion according to the corresponding probability amplitudes of the initial W-type state.This gives another n s − 1 conditions where the index j 1 and j 2 can be arbitrary chosen from 1 to n s .
Taking n s = 2 as an example, the initial two-mode entangled state (W-type state) takes the form as |ψ s = C 1 |10 +C 2 |01 .So the two restrictions on the relative coupling strengths in the sender side read The relation between k 1 and k 2 thus can be derived Next, by substituting the expression above of k 1 into Eq.( 19), we can get the equation for k 2 There must exist a positive solution for k 2 2 .So we can determine the value for the coupling strengths in the sender side to achieve a successful transfer of multi-mode Wtype states with high fidelity.For more general multi-mode states, the possible existence of multi-excitations makes QST a doable task in principle, but a harder challenge in practice.One possible solution is to go beyond the framework of time-independent rectangle pulse control, and combine our complete dynamic expression in Appendix A.2 with the time-dependent pulse optimization method [57,58].

Entanglement Preparation
Quantum entanglement serves as a fundamental element for quantum computation and quantum communication [59][60][61].The generation and distribution of entangled states have long been the main topic in the field of quantum information processing [12,14,[38][39][40][41].In linear optical systems, entanglement can be prepared by certain bosonic transformation such as beam splitter.Making use of the result of linear bosonic transformation obtained before, we raise a new approach to generate multi-mode entangled states in strongly coupled bosonic systems, and give the shortest possible generation time under our protocol.
We first consider the simplest three-mode model with two boson modes a 1,2 coupled to a single-mode channel c.According to the dynamic solutions Eq. ( 3), the intermediate mode c at an arbitrary evolution time t reads where k 1,2 denote the coupling weights between modes a 1,2 and mode c.Similarly, in order to avoid the error introduced by strong coupling, the condition of conserved total excitations at discrete times is assumed to be satisfied Here, m = 2, 3, . . . is an arbitrary positive integer.We can solve the above constraints analytically and give the similar pulse condition for EP tasks, The time τ and the effective coupling strength g ′ can be easily obtained from Eq. ( 25) as θ ≡ ωτ and ζ ≡ ω/g ′ .By substituting these solutions into Eq.( 23), we can simplify the expression of mode c, Notice that the pulse ansatz Eq. ( 25) is similar to that of QST Eq. ( 7) by scaling the parameter m → m/2.In fact, if one sets an optimized pulse amplitude for EP, the amplitude of mode c(0) is distributed into modes a 1 (τ ) and a 2 (τ ), and the state of mode a 1,2 (0) is transferred to the intermediate mode c(τ ) after an optimized duration of EP.By doing so, an EP task is accomplished.If one fixes the coupling strength and extends the time by another EP pulse duration, the amplitudes in modes a 1 (τ ) and a 2 (τ ) will return back to mode c(2τ ) gradually, and the initial state of modes a 1,2 (0) stored in the intermediate mode c(τ ) will transfer to a 2,1 (2τ ).Then a QST task is completed.In addition, similar to the QST cases, the creation and annihilation operators are also decoupled under this condition, since the compact linear transformation Ψ(τ ) = U A (τ )Ψ(0) is also valid with Assuming that we encode a Fock state |1 on mode c, and leave modes a 1 and a 2 in the vacuum states |0 , it is easily checked that the final state of modes a 1 and a 2 after one pulse of Eq. ( 25) is . Thus, the states with different degrees of entanglement can be obtained by varying the relative coupling weights k 1 /k 2 .Specially, if the coupling weights are set as k 1 = ±k 2 , a Bell state |ψ Bell = (|1 1 0 2 ± |0 1 1 2 )/ √ 2 can be generated in a fast speed between modes a 1 and a 2 in the strong-coupling regime.Then, we introduce the logarithmic negativity E N (ρ) ≡ log 2 ( ρ T 1 ) to quantify the degree of entanglement of the final state [62], where the superscript T means partial transportation of the 2-mode density matrix and ||•|| 1 stands for the trace norm of matrix.The maximum of logarithmic negativity for a 2-qubit system is E N (ρ) = 1, which is achieved for Bell states.In Fig. 4(a) we display the evolution of logarithmic negativity of modes a 1 and a 2 for m = 2, 4, and 6, with equal coupling strengths k 1 = k 2 .It is shown that a maximally entangled state will be prepared after applying the pulse introduced in Eq. (25).Besides, we also define the fidelity between the final state ρ f and the Bell state |ψ Bell , f ≡ tr[|ψ Bell ψ Bell | tr c (ρ f )], and show the results of 1 − f for m = 2, 3, . . ., 7 in Fig. 4(b).It is obvious that with our optimized pulse ansatz, the EP fidelity is always equal to 1 for arbitrary cases with m > 1.For comparison, the RWA method always gives a non-negligible infidelity, which can reach as high as 1% in the strong-coupling regime.Thus, our optimized method provide a promising EP scheme to prepare perfect entangled states with low time cost.
Remarkably, considering that in the optimized scheme the EP time is monotonically decreased with increasing coupling strength (and thus a smaller m), we conclude that the fastest time possible to generate a perfect entanglement can be obtained by setting m = 2, which gives ζ = ω/g ′ = 10  3 , an θ = ωτ = √ 10 2 π.An even smaller choice of m < 2 corresponds to an unbounded potential trap and the sine functions with imaginary variant in Eq. ( 23) are replaced by hyperbolic sine functions, leading to a diverging excitation number and an inevitable error (see Appendix B for details).
Next, we show the above method can be directly generalized to the preparation of N-mode entangled states from a 1,2,...,N , where the evolution of the intermediate mode c can be expressed as By defining the united mode a(t) ≡ n j=1 k j a j (t) √ n j=1 k 2 j as in Appendix A.2, the multi-mode EP scheme can be viewed as a linear transformation between c(t) and the boson modes a(t)'s, expressed as Thus, the N-mode entangled state can be obtained as where |0j ≡ ⊗ i =j |0 i .Specifically, if the amplitudes k j are set equal, an N-mode Wtype state [42] In addition, we notice that the proposed N-mode linear bosonic transformation for EP shown in Eq. ( 29) takes the same form as a two-mode EP process except a different definition of a(t).Thus, we conclude that the minimum EP time is independent on the number N of modes.

Summary
We introduce a universal optimized strategy for realizing fast linear bosonic transformation in the strong-coupling regime, where the creation and annihilation subspaces are decoupled, thus preserving the total excitation number and suppressing the infidelity of various quantum information tasks.Based on this strategy, we obtain the optimized pulse ansatz for quantum state transfer (QST) and entanglement preparation (EP), and demonstrate that both tasks can be achieved with high fidelity and fast speed.
Firstly, we simulate the QST task between two boson modes coupled through a single-mode intermediate channel with rectangle shaped pulses, and demonstrate a reduction of infidelity up to 20% in the strong-coupling regime.By analytically solving the system, we obtain a tradeoff relation between the transferring speed and the tolerable error.This result can facilitate the choice of pulse amplitude and duration in experiment to optimally balance the QST time and fidelity.The proposed scheme is completely immune to thermal noise of the channel, and is robust against fluctuations of pulses.It can be further improved to approach higher fidelity by applying a local rotation of the final state to compensate the phase shift induced by the strong coupling, which is also universal and can be obtained analytically.Then we generalize the proposed method to QST of an arbitrary multi-mode W-type state through the common channel mode, demonstrating the potential to extend our protocol to more general quantum information tasks involving multi-mode dynamic processes.Secondly, we apply the optimized method to EP and propose a new approach to prepare multi-mode entangled states including Bell states and W-type states.The fastest preparation time based on our scheme has been obtained, which is independent on the number of the entangled modes.Our results provide new possibilities to achieve high fidelity and fast speed in QST and EP with realistic experimental techniques, and help reduce the stringent requirement of coherence time and temperature that would otherwise be required in quantum computation and quantum communication.In addition, the method may be combined with pulse optimization approach [57,58] to give even better pulse design.
With that, the equations of motion can be rewritten as the following ordinary differential equations (ODEs), This set of first-order ODEs can be solved via a general approach as described below.The vector − → Y (0).In particular, the dynamic solution of the i-th mode can be achieved as follows (Einstein convention is applied here), Here we rearrange the i-th row of the transformation matrix S −1 into the diagonal form diag(S −1 i ), such that it can interchange position with the oscillation matrix e Λt .Using the general procedure outlined above, the four eigenvalues of the ODEs (A.5) are ±i ω 2 − 2gω and ±i ω 2 + 2gω.According to the initial conditions, we get the exact solution of the dynamic evolution of mode a, where ζ = ω/g.The solution for mode c can be obtained analogously, just by exchanging a and c in the solution of mode a.Notice that in the case of |g| > ω/2, we define to avoid possible confusion of double-valued square root.Some special attention needs to be paid for the case of g = ±ω/2, where the denominators of several elements in the matrix approach zero.To overcome this technical difficulty, we remind that the mapping from coupling strength g to the dynamic evolution function should be continuous and analytic.Thus, we can take limitation to calculate the expression Eq. (A.7) for g → ±ω/2.For example, the result for the special point g → ω/2 reads For a small coupling strength g, or equivalently a large value of ζ, the expression above can be rewritten as Retain the first-order term of the frequency ω 2 ± 2gω ≈ ω ± g, and set gt = π/2, we can get the identical results with RWA,

. Multiple modes coupled to an intermediate channel
The more general case with multiple modes coupled to a single intermediate channel mode can be solved analogously.The Hamiltonian reads By defining the united mode a(t) , and rearranging its distribution over individual boson modes, one can obtain the solutions of dynamic evolution as Fast quantum state transfer and entanglement preparation in strongly coupled bosonic systems19 where ζ = ω/g ′ and g ′ = n j=1 k 2 j g.Therefore, one can predict arbitrary correlation functions exactly.

Appendix B. Optimized scheme of quantum state transfer
The conventional quantum state transfer (QST) scenario obtained under RWA suggests the pulse condition g ′ τ = π, under which the fidelity of the task drops significantly in the strong-coupling regime.The exact solution obtained in the previous section can provide an optimized condition with better performance.
For the case of two bosonic nodes coupled to an intermediate channel with the same coupling strength, the exact solution of mode a 1 is expressed without any approximation as, In the weak-coupling limit where RWA is valid (ζ → +∞, ω 2 ± 2g ′ ω ≈ ω ± g ′ ), the solution predicts the final mode a 1 (τ ) = −a 2 (0), i.e., the QST is ideal.However, in the strong-coupling regime where RWA breaks down, the final mode a 1 (τ ) deviates from a 2 (0) and causes an accumulation of infidelity.According to Eq. (B.1), we find the optimized pulse ansatz to achieve better QST performance as follows, where θ = ωτ is a dimensionless parameter dependent on the system frequency ω and the QST time τ .The parameter ζ = ω/g ′ determines the relative coupling strength.The pulse parameter m is an integer which needs to be chosen appropriately in experiments.Under the optimized pulse condition above, we can easily check that the final mode can be simplified as The parameters ζ and θ can be solved from Eq. (B.2) as This solution gives the optimized pulse scheme of QST.It is easy to find that a larger value of m corresponds to a weaker coupling strength g ′ and a longer QST time τ .Specifically, a strong-coupling regime with g ′ 0.1ω is reached for m 11.
To illustrate the model more explicitly, we write down the potential energy part of the Hamiltonian with canonical variables as . This potential can be expressed in a diagonal form using eigenvectors {Y 1 , Y 2 , Y 3 }, leading to V = i ωi Y 2 i with ω1,2,3 = ω, ω ± 2g ′ .By taking the kinetic energy part T = ω 2 (P 2 1 + P 2 2 + P 2 c ) into account, the Hamiltonian can be diagonalized with the eigenmodes as whose eigenfrequencies √ ω ωi are embodied in the dynamic solution (B.1).In the regime of 2 |g ′ | < ω, i.e., m > 2, the potential trap forms a three-dimensional parabolic surface opening up in all three directions, which establishes a bounded system with positive eigenfrequencies √ ω ωi .While entering into the regime of 2 |g ′ | ≤ ω, i.e., m = 1 or 2, one dimension of the parabolic surface turns to curve down.Such unbounded potential will cause unphysical imaginary eigenfrequency, which we should avoid.

Appendix C. Conservation of the total excitations
We notice that using the optimized scheme of pulse, the operators {a 1 (τ ), a 2 (τ ), c(τ )} can be mapped into the subspace spanned only by the initial annihilation operators {a 1 (0), a 2 (0), c(0)} (without the creation operators).Thus, we can define a transfer matrix S to represent the map from the initial modes to the final modes a i (τ ) =   We first consider a task of transferring a Fock state |n 1 from mode a 1 to a 2 via an ideal channel at zero temperature.After the QST process is completed at time τ , the received state can be expressed as,

21
). (S25) The evolution operator is represented as U (τ ) = T exp[−i τ 0 dt Hint (t)], where T stands for the time-ordered operator and the Hamiltonian is given by Eq. (S12).
To write down the second line of the expression above, we notice that using the optimized scheme of pulse, the operators {a 1 (τ ), a 2 (τ ), c(τ )} can be mapped into the subspace spanned only by the initial annihilation operators {a 1 (0), a 2 (0), c(0)} (without the creation operators).Thus, we can define a transfer matrix S to represent the map from initial modes to the final modes a i (τ ) = j S ij a j (0).Apparently, the S matrix is unitary according to the Using this notation, we can write down two identities j S ij a j (0).Apparently, the S matrix is unitary according to the commutation relations Then, we can easily check the conservation of total excitations by calculating the expectation value Therefore, the carefully selected pulse ansatz makes the creation and annihilation operator effectively decoupled, and the total excitation number can be preserved regardless of the breakdown of the U(1) symmetry.The change of the total excitations for cases of m = 5, 8, and 16 are demonstrated in Fig. C1(a), C1(b) and C1(c), respectively, where we define Ntot ≡ a † 1 a 1 + a † 2 a 2 + c † c.The blue solid (orange dashed) line represents the dynamics from the Hamiltonian with (without) the counterrotation terms.As the coupling strength is weak enough, such as the m = 16 case, the orange dashed line deviates only very little from the blue solid line, indicating that the RWA can be safely applied.However, the total number of excitations manifest greater and non-negligible fluctuations in the strong-coupling regime for m = 5.In another aspect, the fluctuation exhibits obvious periodicity.And the number of excitations keeps unchanged after a complete cycle.We set our QST condition right after one cycle to suppress the error brought by the fluctuation of excitations.
The evolution operator is represented as , where T stands for the time-ordered operator and the Hamiltonian is given by Eq. (1).
Using the notation defined in the previous section, we can write down two identities where S T is the transpose of S. Thus, the final state U(τ ) |0 1 , 0 c , 0 2 can be obtained as

This expression then leads to another relation
Here, |n i denotes the Fock state of mode O i , |ψ n i is the state in the orthogonal subspace spanned by O j =i , and A n i is the corresponding coefficient.Owing to the orthogonality of the Fock basis, we conclude that A n i =0 = 0 for all three modes O i , and thus U(τ In the last line of Eq. (D.1), we expand the expression of the final state in orders of K 11 /K 21 , since the conditions of K 21 ∼ 1 and K 11 ∼ 0 are naturally expected to obtain high fidelity.To the first order, the error is A straightforward estimation of the error in general cases gives 1 Using the definition Eq. (B.1), we obtain the exact result of K 11 as In Eq. ( 11), the tradeoff relation is expressed as the constraint function F (ρ, 1 − f, τ ) = 0. Obviously, the function can be acquired according to the analysis According to the calculation above, the coefficient The addition of the two black vectors gives the coefficient K 21 (green vector), while the difference gives the coefficient K 11 (red vector).As the parameter m gets smaller, the upper black vector will rotate counterclockwise (the direction of the red thick arrow), making the green vector rotate counterclockwise and leading to a larger value of θ r .
2)θ(m)/4], where m is viewed as a continuous parameter temporarily.One can easily check that G(m) is a monotonically decreasing function of m.
Once the error that the receiver can tolerate is set as E tol , the QST time τ can be evaluated by Eq. (D.8), (D.9) Equivalently, it can be rewritten as G(m) < E tol / n 1 ρ .By solving the inverse function of G(m), we can get the threshold of parameter m, which reads m > m th ≡ G −1 E tol / n 1 ρ .Thus, the threshold of τ can be obtained by θ th = θ(⌊m th ⌋ + 1), where ⌊m th ⌋ is the greatest integer less than or equal to m th .Since θ th ≡ ωτ th , the shortest transfer time is given by τ th = θ th /ω.For instance, the results for transferring a Fock state |n 1 with n 1 = 1, 2, 3 via a zero-temperature channel is displayed in Fig. D1(a), while the ones for a coherent state |α with |α| = 0.6 ∼ 2.0 is shown in Fig. D1(b).Here the analytic results are calculated with the expression Eq. (D.8), while the numerical results are obtained via a direct evolution of the equations of motion associated with the Hamiltonian without RWA.Notice that the analytic prediction agrees perfectly with the numerical results in all parameter regimes.

Appendix E. The extra rotation in phase space
Strong coupling effect induces a non-negligible rotation of the initial state in phase space.Using our optimized pulse ansatz Eq. (B.3), we can obtain the final state, ) is the intrinsic result associated with the energy shift from ω ± g ′ in the weak-coupling regime to ω 2 ± 2g ′ ω in the strong-coupling regime.The second term is the error deviated from the desired mode a 2 (0), and the coefficient is given exactly as in Eq. (D.7).We can calculate the coefficients above according to a geometry relation in the complex plane.As illustrated in Fig. D2, the rotation angle θ r increases with the decreasing of QST time θ = ωτ .Such behavior can be understood by the red arrow in the right subfigure.The black vector lying on the real axis stands for the number 1/2, and the other black vector above the real axis refers to (1/2)e −i[( √ 1+2/ζ+ √ 1−2/ζ)/2−1]ωt = (1/2)e 2iθr .The superposition of the two vectors gives the coefficient cos(θ r )e iθr corresponding to the green vector in the subfigure, and the difference of the two vectors (red vector) gives the coefficients K 11 corresponding to the error of the QST task.Apparently, as the parameter m gets smaller, the QST time becomes smaller.As a result, the second black rector will rotate more counterclockwise, making the length of the red vector (error) become lager.Meanwhile, the green vector will rotate more counterclockwise, leading to a larger angular shift in phase space.
With that, an additional local rotation applied to node a 2 is proposed after the optimized QST pulse of Eq. (B.3), in order to further suppress the infidelity caused by the rotation θ r .It can be realized by implementing a local pulse on node a 2 expressed as H r = g r (t)a † 2 a 2 with the pulse area g r (t)dt = θ r .Here, we take an example of transferring a coherent state |ψ = |e iπ/2 in the strong-coupling regime g ′ 0.1ω (corresponds to 5 ≤ m ≤ 11 in our optimized QST method).As shown in Fig. E1, such a local rotation, labeled as "further optimization", can suppress the infidelity for 2 ∼ 3 orders of magnitude, which is quite remarkable considering the error is already less than a few percent.We emphasize that this further optimization method is efficient for coherent states since they are sensitive to the rotated phases in phase space.Thus, it is also useful to suppress the QST error when transferring a general state that is sensitive to phase, such as squeezed states, cat states, etc.On the other hand, this method makes no difference to Fock states, which are centrally symmetric in phase space and thus insensitive to a phase rotation.

FIG. 2 . 2 (
FIG. 2. (a) The infidelity 1 − f versus the channel temperature T with the parameter k = 6, for transferring the |1 (red) and coherent state |e i π 2 (blue) by the two methods respectively.(b) The infidelity for transferring the Foc with different deviations of the pulse time |∆τ | in the optimized pulse ansatz.(c) the infidelity changes with the pha coherent state with fixed amplitude |α| = 1.0,where the pulse parameter is set as k = 5.In (b) and (c), a zero-te channel T = 0 is assumed.

Figure 2 .
Figure 2. (a) The infidelity versus channel temperature T for tasks of transferring the Fock state |1 (red) and the coherent state |e iπ/2 (blue) by the optimized and RWA methods, respectively.The pulse parameter is chosen as m = 6.(b) The infidelity for transferring the Fock state |1 with different pulse duration fluctuations |∆τ | in the optimized pulse scheme.Here, an ideal channel with T = 0 is assumed, and solid lines as guides for the eyes are drawn to connect symbols.

FIG. 3 .
FIG. 3. The initial (top row) and final (bottom row) states of a QST of CAT state |CAT(α) + with α = 1.2 through the optimized method.The pulse parameter is chosen as k = 11.(a) When the local rotation is not applied, the fidelity of the final state is f = 0.9819.(b) After the local rotation, the fidelity is increased to f = 0.9922.

Figure 3 .
Figure 3.The initial (top row) and final (bottom row) states of a QST of cat state |cat(α) + with α = 1.2 through the optimized method.The pulse parameter is chosen as m = 11.(a) When the local rotation is not applied, the fidelity of the final state is f = 0.9819.(b) After the local rotation, the fidelity is increased to f = 0.9922.

Figure S5 .Figure 4 .
Figure S5.The task to generate the maximum entangled state in the nodes 1,2 with the equal coupling strengths k 1 = k 2 .(a) The Change trend of the log negativity for cases k = 2, 4, 6.(b) Comparison between our optimized method with the RWA method in the fidelity of the final entangled state for the cases k = 2, 3, . . ., 7.
as the function of time t is driven by the coefficient matrix K, i.e., − → Y = K − → Y .Diagonalizing the matrix K = S −1 ΛS, we can reorganize the modes into the eigenmodes − → Y ′ with independent oscillation frequencies − → Y ′ = Λ − → Y ′ .Thus, the solution of the reorganized modes can be obtained directly as − → Y ′ (t) = e Λt − → Y ′ (0).The solution of the original modes can then be expressed as − → Y (t) = S −1 e Λt S 0

Figure S1 .
Figure S1.The total number of excitations Ntot changes against the time t for different selections of k, (a) k = 5; (b) k = 11; (c) k = 16.In the bottom half, we arbitrarily select the coupling strengths in the strong coupling regime as (d) g = 0.15ω; (e) g = 0.17ω; (f) g = 0.2ω.Here the blue solid(orange dashed) lines stands for the Hamiltonian with(without) the counterrotation terms.The initial excitations number is set as Ntot = 1 in all cases.

Figure C1 .
Figure C1.The total number of excitations Ntot changes against time t for different selections of m, (a) m = 5; (b) m = 8; (c) m = 16.Here the blue solid (orange dashed) lines stands for the Hamiltonian with (without) the counterrotation terms.The initial excitations number is set as Ntot = 1 in all cases.

Figure D1 .
Figure D1.The tradeoff relation between the evolution time τ and the infidelity 1 − f .(a) Transferring Fock states |1 (red), |2 (green), |3 (blue).The solid lines represent the analytical results of the error, while the dots show numerical results of dynamic evolution using our optimized QST method.(b) Transferring coherent states with amplitude |α| = 0.6 ∼ 2.0.The curved surface shows the analytical relation, while the dots around the surface are the numerical results.The range of m is chosen from 5 ∼ 17.

Figure D2 .
Figure D2.The left figure displays the rotation angle θ r changes against the QST time τ .The right figure illustrates the phase error originated from strong coupling.The addition of the two black vectors gives the coefficient K 21 (green vector), while the difference gives the coefficient K 11 (red vector).As the parameter m gets smaller, the upper black vector will rotate counterclockwise (the direction of the red thick arrow), making the green vector rotate counterclockwise and leading to a larger value of θ r .

Figure E1 .
Figure E1.The infidelity when transferring a coherent state |ψ = |e iπ/2 through an ideal channel in the strong-coupling regime with 5 ≤ m ≤ 11.It is found that the optimized method without further optimization (blue solid line) suppresses the infidelity in the conventional RWA method (orange dashed line), while a more remarkable decrease is observed after applying the further optimization with local rotation (green dashed line).
ψ| being the density matrix of the sending state at node a 1 and ρ f ≡ tr 1,c [ρ 1,c,2 (τ )] the reduced density matrix of the received state at node a 2