Quantum trajectories and their statistics for remotely entangled quantum bits

We experimentally and theoretically investigate the quantum trajectories of jointly monitored transmon qubits embedded in spatially separated microwave cavities. Using nearly quantum-noise limited superconducting amplifiers and an optimized setup to reduce signal loss between cavities, we can efficiently track measurement-induced entanglement generation as a continuous process for single realizations of the experiment. The quantum trajectories of transmon qubits naturally split into low and high entanglement classes corresponding to half-parity collapse. The distribution of concurrence is found at any given time and we explore the dynamics of entanglement creation in the state space. The distribution exhibits a sharp cut-off in the high concurrence limit, defining a maximal concurrence boundary. The most likely paths of the qubits' trajectories are also investigated, resulting in three probable paths, gradually projecting the system to two even subspaces and an odd subspace. We also investigate the most likely time for the individual trajectories to reach their most entangled state, and find that there are two solutions for the local maximum, corresponding to the low and high entanglement routes. The theoretical predictions show excellent agreement with the experimental entangled qubit trajectory data.


I. INTRODUCTION
Measurement-induced entanglement of spatially separated quantum systems is a startling prediction of quantum mechanics [1][2][3][4][5][6][7]. Recent experiments have demonstrated this effect via single-photon heralding [8][9][10][11], as well as via continuous measurement of photons interacting with qubits [12]. The latter has the advantage of being able to investigate the physics of entanglement creation continuously, leading to new effects such as the sudden creation of entanglement after a finite measurement period, dubbed entanglement genesis [6]. However, many questions are outstanding, such as (a) What is the complete characterization of the dynamics of entanglement creation as a continuous trajectory?, (b) What is the statistical distribution of the entanglement at any time during the process?, and (c) What is the most likely way entanglement is generated? In this work, we give a systematic answer to these questions and more, by analyzing experimentally entangled quantum trajectories of jointly measured transmon qubits and showing excellent agreement with the theory developed here.
The rapid development of quantum information science in the superconducting domain [13] has seen an exponential increase in qubit coherence time within the past decade, leading to many scientific advances [14]. This technological progress has led to a wide variety of advances in quantum physics as observed and controlled in these systems, including > 99% fidelity in single qubit quantum gates [15], multi-qubit entanglement generating gates [16], the violation of Bell's inequality [17] and quantum process tomography [18]. Recent developments include the observation of quantum states of light in resonators [19], as well as nearly quantum-limited parametric amplifiers [20,21].
The improvement of coherent quantum hardware has brought with it a renewed focus on the physics of quantum measurement. Generalized measurements have been carried out in superconducting qubits [22], realizing probabilistic measurement reversal [23,24], weak values [25,26], and their connection with generalized Leggett-Garg inequalities [27]. Continuous measurements [28] in superconducting systems have only recently been realized, owing to the challenge associated with high-fidelity detection of microwave signals near the single-photon level. In particular, experimental achievements include continuous feedback control [29,30], the tracking of trajectories in individual experiments in both the plain measurement case [31][32][33][34], as well as with a concurrent Rabi drive [35]. These experiments show detailed quantitative agreement with theory, indicating good understanding of quantities such as the most likely path of the quantum state between boundary conditions, predicted with an action principle of an associated stochastic path integral [36][37][38].
Going beyond a single qubit opens the possibility of measurement-induced entanglement, using a dissipative process as a tool to generate quantum correlations. For quantum architectures building upon transitions at optical frequencies [8][9][10][11], this feat has typically been achieved by relying on the correlated detection of photons at the output of a beam-splitter [1] to herald an entangled state. While powerful, this measurement protocol is binary and instantaneous, and allows no insight into the dynamical processes underlying the generation of the entangled state. In solid state systems, such as superconducting qubit transitions in the microwave regime, there has been tremendous interest in continuously generating bipartite [3][4][5][6] and multi-partite [2,39,40] entangled states, using weak measurements that slowly interact with the qubits, in such a way that enables the resolution of the dynamical aspects of the entangling backaction. Analog feedback control can be naturally applied in the weak continuous measurement regimes [41][42][43] and digital feedback-generated entanglement has already been demonstrated [44]. Joint measurement is uniquely useful as a means to generate entanglement between remote qubits [12,[45][46][47][48][49], for which no local coupling exists and therefore no unitary means of generating entanglement are available.
The chief advantage of the continuous approach is in the efficiency of entanglement generation. In contrast to photon-counting schemes (in which entanglement generation rates are heavily limited by photon losses), continuous measurement enables an entanglement generation rate limited only by the pre-measurement initial state and the postselection criterion, which are both experimental choices. However, continuous measurements may be highly sensitive to dissipative losses and inefficiencies in amplification. Understanding and studying the dynamical processes underlying continuous measurementinduced entanglement is therefore critical for balancing these tradeoffs.
In this work, we combine an efficient quantum amplifier with a continuous half-parity measurement, to conduct detailed experimental and theoretical investigations of the statistics of individual trajectories of qubit pairs as they undergo the entangling process. By peering into the ensemble, we can understand the full spectrum of evolution paths as the two-qubit state gradually projects onto the entangled subspace or onto a trivial separable state. We explore the probability distribution of the qubits' concurrence to understand how the distribution changes in time, from a separable state with zero concurrence, to projected states in either a separable subspace or an entangled subspace. This can also be seen from the most likely path analysis, showing the emergence of different most probable paths for each final state. Moreover, we investigate the distribution of the time-tomaximum-concurrence, finding that the most probable time to maximum values of concurrence has a bimodal structure. Studying the statistics of a large set of trajectories -rather than averaging over all of them to study the dynamics of the ensemble -enables an unprecedented understanding of the dynamics of entanglement creation under measurement.
The paper is organized as follows. In Section II, we describe our transmon qubit measurement setup and the method of reconstructing the joint trajectories. In Section III, we derive the concurrence-readout relation and compute the distribution of concurrence as the measurement back-action proceeds. This distribution is then compared with the data generated from the experiment. In Section IV, we turn to the most likely path analysis, finding three possible paths that the joint system takes from a separable initial state to final entangled or separable states. The experimental most likely paths are generated independently to compare with the theoretical prediction. We also discuss the distribution of the time for the two qubits to reach their maximum entanglement. The conclusions are presented in Section V. Additional details of the most likely path calculations and a discussion of the parity meter are presented in the Appendices.

II. TRAJECTORIES OF TRANSMON QUBITS
We consider a two-qubit system realized by superconducting transmon qubits embedded in spatially separated microwave cavities in a setup optimized to reduce losses between the cavities. These qubits are jointly measured via a dispersive readout in a bounce-bounce geometry (see figure 1(a)) in which a microwave tone is sequentially reflected off of two copper cavities each containing a transmon qubit, and subsequently is amplified and measured via homodyne detection. The cavities are directly joined by a circulator that enforces the unidirectional transfer of the coherent state. A single transmon gives a dispersive phase shift of ±φ for states 0⟩, 1⟩. The shift φ is given in general by φ = arg α 0⟩ − arg α 1⟩ , where α i represents the intra-cavity coherent state conditioned on qubit state i. In the fully symmetric, weak measurement case (when all qubit and cavity parameters are identical and χ ≪ κ), this phase shift is given by φ ≈ 2χ κ, where χ is the cross-Kerr nonlinearity between the qubit and the cavity, and κ is the cavity damping rate.
In the bounce-bounce geometry, there generally exists a probe frequency at which the dispersive shift from both transmon qubits are the same (φ 1 = φ 2 = ϕ), such that the measurement tone can acquire a phase shift of either 2ϕ, 0, 0, −2ϕ for states 00⟩, 01⟩, 10⟩, 11⟩. For small χ κ, this results in a half-parity measurement on the two qubits, where the measurement result can distinguish between three subspaces: the 00⟩ state, the 11⟩ state, and the odd-parity subspace, but not within the odd parity subspace, spanned by the 01⟩ and 10⟩ states. Conditioning on the measurement results with a zero phase shift can lead to creation of an entangled state within the odd parity subspace (superpositions of 01⟩ and 10⟩).
When the output of the second cavity is directed to a nearly-quantum limited amplification chain, the instantaneous homodyne detection signal can be correlated with the measurement back-action on the qubits, and therefore can be used to track the evolution of the system in time. By reducing the amplitude of the coherent state used to measure the system, we can engineer an entangling measurement with a characteristic measurement time ranging from several hundred nanoseconds to several microseconds: critically, these timescales are easily resolvable experimentally. The dynamics or the trajectory

Matrix elements
Matrix elements x 1 x 2 x 3 x 4 x 5 x 1 x 2 x 3 x 4 x 5 of the system state can be obtained via the full master equations [45,46], using a two-cavity polaron transformation to account for the cavity degree of freedom, giving the stochastic master equation for the qubit trajectories. Alternatively, in a limit of large cavity decay rate κ ≫ χ , the qubits evolution can be continuously tracked via the quantum Bayesian approach [2,3], inferring the current states of the system from the measurement readouts and how likely they are to occur. The two approaches both show good agreement in tracking the qubit pair state [12].
In this paper, we focus on the quantum Bayesian approach, as it is directly related to the probability distribution of the measurement readout, and naturally leads to the probability distribution of quantum trajectories. Let us denote p(V t i) as a probability density function of a measurement readout V t conditioned on the twoqubit states i, where i = 1, 2, 3, 4 representing the states 00⟩, 01⟩, 10⟩, 11⟩ respectively. The quantum Bayesian update for this type of double-qubit measurement provides a convenient way to calculate the joint state at time t, given a known state at the initial time and the readout V t , where ρ ij denotes the ij element of the two-qubit density matrix and γ ij is a decoherence rate associated to the matrix element. We define the readout V t ≡ (f t) ∫ t 0Ṽ (t ′ )dt ′ − v 0 as a time-average of a raw homodyne voltage signal rescaled with a weight factor f and an off-set v 0 , where f is chosen so that the variance σ 2 Vt = 1 4η m t is a function of a quantum efficiency of the homodyne measurement, η m ≈ 0.22. The total probability distribution p(V t ) = ∑ 4 k=1 ρ kk (0)p(V t k), shown in Figure 1(b), slowly resolves into the three peaks expected for a half-parity measurement.
The conditional readout distributions are wellapproximated by Gaussian functions, giving p( The measurement process cannot distinguish the two states in the odd-parity subspace, therefore the readout distributions corresponding to the states 01⟩ and 10⟩ are completely (or nearly) overlapped, giving δv 2 ≈ δv 3 ≈ 0 and −δv 1 ≈ δv 4 ≈ δv. The measurement strength is characterized by an inverse of a characteristic measurement time τ m ≈ 1 δv 2 η m . The dephasing rates γ ij for δv i ≠ δv j are dominated by the effect of the distinguishability between states i and j, γ ij ∼ (η −1 m −1)(δv i −δv j ) 2 4s [2], resulting in the strong suppression of all off-diagonal elements except ρ 23 . In an ideal half-parity measurement, the decay of ρ 23 would be limited only by the intrinsic lifetimes of the qubits; however, we must additionally account for experimental imperfections in the matching of δv 2 and δv 3 and for the loss of photons between the two cavities. These effects are included in the (slightly time-dependent) dephasing rate γ 23 .
Since we expect most of the off-diagonal terms to damp quickly, we only consider five density matrix elements: 44 , and x 5 ≡ ρ 23 [50]. In order to compare the Bayesian prediction to the true density matrix, we perform conditional tomography [12,51] in order to experimentally reconstruct the full experimental mapping V t ↦ ρ(V t , t). In Figure 1(c), we show good agreement between the Bayesian prediction and conditional tomography reconstruction of x i (V t ) for a single measurement time t = 0.48 µs. We show examples of transmon qubit trajectories in Figure 1(d) for an initial state prepared in a product ofx-states such that We show both the Bayesian reconstruction and the tomographic verification, which are in good agreement. This verifies that the Bayesian reconstruction can be used to faithfully translate a noisy measurement signal into the stochastic evolution of a joint quantum state.

III. CONCURRENCE TRAJECTORIES AND THEIR DISTRIBUTION
As a measure of entanglement between two parties such as the transmon qubits, concurrence is a convenient choice and can be computed directly from the density matrix of the system [52]. The concurrence formula for the half-parity setup is greatly simplified because of the suppression of most matrix elements, resulting in a X-shape density matrix, of which the concurrence is calculated from [52,53], The value of concurrence ranges from 0 for a separable state to 1 for the Bell states. The concurrence trajectory of an exemplar trajectory is shown in the inset to Figure 1(d). In this section, we will use the simplified formula Eq. (2) to show that the concurrence of the twoqubit state can be determined directly from the measurement readout, which then leads to the derivation of the concurrence probability distribution as a function of the measurement readout and measuring time.

A. Concurrence-readout relationship
Let us consider the concurrence formula in Eq. (2), where its value is determined by the second term in the bracket, which we denote c t ≡ 2 If c t is a non-negative quantity (i.e., c t ≥ 0), then the concurrence is simply given by C(t) = c t . We will show at the end of this section that this is always the case for our chosen initial qubit state and parameter regimes, but is not true in general [6]. From the Bayesian update in Eq. (1), we calculate the quantity, where x 0 i for i = 1, ..., 5 are the matrix elements of the initial qubits state, and γ = γ 23 . We have used simplified notations P i ≡ p(V t i) for the probability distributions, and N for a normalized factor given by N = ∑ 4 k=1 x k (0)P k . Substituting the probability distribution functions with the Gaussian functions of the means δv i for i = 1, ..., 4, we obtain a form of c t explicitly as a function of V t and t, where the prefactor is given by The quantity c t in Eq. (4) would represent the actual concurrence of the qubits state at any time t, if c(t) ≥ 0 is satisfied. For our chosen initial state, a product of single qubitx-states: From the experimental data (e.g., for the setup with τ m = 0.60µs), we have −δv 1 ≈ δv 4 and δv 2 ≈ δv 3 ≈ 0 (giving α 14 , α 23 , β 23 ≈ 0) and β 14 ∼ 3.2 MHz, while γ < 0.6 MHz. Therefore, the condition is always satisfied and the second term in the bracket of Eq. (4) decays faster than the first term, resulting in an always non-negative quantity. Consequently, the quantity in Eq. (4) gives the concurrence-readout relationship C(V t , t) = c t (V t , t), and the concurrence at any time t can be determined directly from the time-averaged measurement readout V t .
The concurrence formula in Eq. (4) can be simplified further by considering a perfectly symmetric half-parity measurement, δv 2 = δv 3 = 0 and −δv 1 = δv 4 = δv. Given the initial state, a product of two qubitx-states, the concurrence is then given by, where the subscript 'ps,x' indicates the perfectly symmetric half-parity measurement given the specific initial state.

B. Probability density function for concurrence trajectories
From the direct relationship between the measurement readout and the concurrence of the qubits state, the probability density function of the concurrence can be derived from the probability distribution of the time-averaged signal V t . The distribution of the time-averaged readout is given by, The variance of the distribution σ 2 Vt = s 2t narrows as time increases, leading to the collapse of the joint qubit state into three categories: 00⟩ state, 11⟩ state, and some superposition state of 01⟩ and 10⟩ after a few characteristic measurement times τ m .
Knowing the probability density function of the timeaveraged signal, we follow the transformation of random variables V t → C using the concurrence-readout relationship Eq. (4) (or (5) for a perfectly symmetric case). The concurrence C(V t , t) is not a monotonic function in V t ; instead it has a bell-like shape as shown in the inset of Figure 2(a). We write the cumulative distribution func- is a probability density function for the concurrence, and V + , V − are two solutions that arise from solving Eq. (4) (or (5) The concurrence distribution is then obtained by taking a derivative of the cumulative distribution, noting that V − (c, t) and V + (c, t) are functions of the concurrence c and time t. The full solution of p C,t (c) is quite lengthy and is not shown. We show in Figure 2(a) the plots of concurrence probability distributions (7) for different values of time t, and in Figure 2(b,c) the density plot comparing between the theory and the transmons experiment. At an early time, the distribution of concurrence is narrowly peaked near its maximum which increases over time, whereas at later times, a second peak emerges near the zero concurrence, showing a bimodal distribution. In Figure 2(b), the theoretical histogram for the concurrence is obtained by integrating the theory curves Eq. (7) for the probability over small intervals δc ≈ 0.015. This is to make a fair comparison with the histogram of the experimental data in Figure 2(c), calculated with a bin size of 0.015. We note that a short delay in the experimental entanglement creation is a result of the cavity ring-up time, which will be discussed more in the next section.
We stress that the concurrence distribution has a sharp upper bound (shown as a grey dotted curve in Figure 2(a)), which the concurrence cannot exceed. In order to understand why the probability distribution for the concurrence has a sharp upper cut-off at any time, we recall that the density matrix of the two-qubit system, conditioned on the time-integrated readout V t , is entirely specified by that (random) outcome, together with the initial state, the dephasing rate, and other parameters of the problem, Eq. (1). Consequently, the concurrence is controlled by V t , as in Eq. (5). As can be seen from the inset of Figure 2(a), the concurrence, plotted as a function of the measured signal V t is bounded from above for any fixed time t by some amount we call C max , and consequently, any value of concurrence higher than that maximum (whose value will change as the time increases) cannot be realized. Therefore the probability distribution of concurrence has a sharp upper cut-off given by C max (t). Physically, this indicates that there is an upper limit on how fast entanglement can be created by the continuous measurement in this situation, even for rare events of the measurement process.
For the perfectly symmetric case in Eq. (5), we can find an analytic solution for the upper bound of the concurrence, knowing that cosh(x) has its minimum at x = 0, C ps,x (V t ) has its maximum at V t = 0, so consequently, the concurrence upper bound is given by The behaviour of this bound is a result of two competing rates, between the extra dephasing rate γ and a measurement rate δv 2 s. Eq. (8) increases from zero for small time and decays for long time after reaching its maximum concurrence as seen in Figure 2. The maximum possible concurrence for this qubit half-parity measurement and the time this happens can be obtained from this relation. More about the time to reach maximum concurrence will be discussed in Section IV C.

IV. MOST LIKELY PATH ANALYSIS
In addition to the distribution analysis in the previous section, where we treated the concurrence at each point in time independent from any other times, we now incorporate the notion of connected trajectories, finding a probability density function for quantum trajectories and their most likely paths. These most likely paths describe routes with highest probability density, taking into account that each of the points in the trajectory ensemble are connected via quantum state update rules, e.g. in Eq. (1). As we have seen previously, the concurrence distribution exhibits a transition from a single-peak distribution to double-peak distributions, one at the upper bound concurrence and another near zero concurrence. Here, with the notion of trajectories, we will also see that the two-qubit state starting in its initial state gradually collapses to three subspaces: the 00⟩ subspace, the 11⟩ subspace and the odd subspace, as described by three most likely routes.
Let us consider the five non-trivial elements of the two-qubit density matrix {x 1 , x 2 , x 3 , x 4 , x 5 }, and treat each element as an independent variable. Each state trajectory is represented by a time series of these five elements, which corresponds to one realization of the measurement readout {v t } = {v 0 , v δt , v 2δt , ..., v t } where we define, v t = (f δt) ∫ t+δt tṼ (t ′ )dt ′ − v 0 , as an instantaneous readout at time t with an integration time δt. Since the readout v t 's are assumed Markovian and only depend on the qubit states right before the measurement, a joint probability density function for the readout realization is given by, a product of probability distributions of v t ′ from t ′ = 0 to t ′ = t with a time step δt. The probability function for an instantaneous readout is given by p(v t i) = δt πs exp{−(v t − δv i ) 2 δt s} for i = 1, 2, 3, 4.
To derive the most likely path for these two-qubit trajectories, we use the Bayesian update equation for the two-qubit state Eq. (1), adapted to a state update every time step δt, and then we maximize the joint probability density Eq. (9) constraining the state update equations. Following the most likely path analysis for quantum states under continuous measurement [36,37] and introducing Lagrange multipliers {p 1 , p 2 , p 3 , p 4 , p 5 } for the contraints, we obtain differential equations for an optimal path in the qubits state space, where i = 1, ..., 4 for the first line and the variables x k are time-dependent functions. We note that v t in these equations behaves as a "smooth" optimal readout, which is a function of the qubit density matrices and their Lagrange multipliers, determining an optimal path from an initial state to its final state [36]. An optimal path is a solution of 10 differential equations: 5 for the qubit state variables Eq. (10), another 5 for Lagrange multipliers, and the optimal readout as a function of both sets of variables (see Appendix A). These equations describe most likely paths for a measurement with any values of δv 1,...,4 and can also be generalized to include the effect of external drives on the two-qubit. However, in the absence of external drive, these can be simplified as we will see in the next section.

A. Most likely paths for joint measurement of transmon qubits
In this work, where the transmon qubits only evolve with the influence of the measurement back-action, the optimal readout is found to be constant in time (see Appendix A). We can therefore bypass solving the full set of the differential equations [48], and only compute the qubit evolution in Eq. (10) with constant v t , looking for measurement results with maximum likelihood. To estimate the likelihood of a two-qubit trajectory, we evaluate a logarithm of the probability density function Eq. (9) approximated to first order in δt [36] to obtain, where S 0 represents a state-independent part of the joint probability function P. We show in the insets of    ure 3(a,b,c) examples of the approximated log-likelihood as functions of optimal readout v t , for three different measurement strengths. In the case of τ m ≈ 2.10µs, there is only one maximum likelihood value located at v t ≈ δv 2,3 , which gives the most likely path with high concurrence shown as a solid curve in Figure 3(a); whereas, in the stronger measurement cases, Figure 3(b) and (c), the approximated log-likelihood have three local maxima: the middle ones, v t ≈ δv 2,3 , corresponding to most likely paths collapsed to entangled states (high-concurrence branches), and the sided peaks, v t ≈ δv 1 and v t ≈ δv 4 , corresponding to two branches of most likely paths collapsed to 00⟩ and 11⟩ states (low-concurrence branches), respectively. We note that this technique of finding multiple most likely paths is not under final-state constraints as in [35,36]. In order to test the theoretical most likely paths, we need to extract the most likely paths from the experimental data. We collect an ensemble of 10 4 transmon trajectories and then compute average trace distances between any two trajectories, (e.g., ρ a and ρ b ), for all possible pairs. The goal is to pick the first few trajectories with minimum total distance to other trajectories in the ensemble, and average them to get an estimate of the experimental most likely paths. For this particular set of data, we choose ∼ 10 2 highly likely trajectories to get a smooth estimate of the most likely paths. However, for the case that there exist multiple (e.g., three) most likely paths, we first divide the ensemble into subensembles according to their trace distance, and then apply the minimum total distance procedure to the trajectories in each subensemble separately. In Figure 3(a,b,c), we show the concurrence of the experimental most likely paths as data points.
Ideally, the theoretical most likely paths predicted from the log-likelihood Eq. (11) for each set of measurement strength using the initial state {x 0 1 , ..., x 0 5 } = {1 4, ..., 1 4} would be good enough to compare with the experimental data. However, in the experiment, after the initial state has been prepared, the cavities take some time to reach their steady state condition, making the parameters δv 1 , ..., δv 4 unstable during the first ∼ 0.13µs time. Therefore, we need to let the initial qubit state evolve and then find new "initial" states at time t = 0.13µs for the theoretical most likely path calculation. We use the experimental most likely states at t = 0.13µs (one for each branch) as the initial states in calculating the log-likelihood functions (the insets), which then lead to an excellent prediction of the most likely paths and their concurrence for the rest of the evolution.
We also note that the unequal two low-concurrence branches in Figure 3(b,c) happen because the population of the states drifts more toward the ground state 00⟩ during the transient time, as a result of the qubit relaxation during the measurement. Moreover, in Figure 3(d,e,f), we present the evolution of matrix elements of the highconcurrence most likely paths, showing a good agreement between the theoretical most likely paths and experimental most likely paths post-selected with the most entangled state at the final time T = 1.6µs.

B. Most likely paths for perfectly symmetric half-parity measurement
We are also interested in finding analytic solutions for the most likely paths for the symmetric case: δv 2 = δv 3 = 0 and −δv 1 = δv 4 = δv. The differential equations (10) for the qubit state simplify to, where we have defined new variables x p ≡ x 1 + x 4 and x m ≡ x 1 − x 4 , and constant parameters a = v t δv s and b = δv 2 4s. The first two equations can be solved independently from the rest. For the case when v t = δv 2,3 = 0, which corresponds to the most likely odd-parity result, we obtain an analytic solution for the high-concurrence branch of the most likely path, where the proportionality factor is the inverse of a normalized factor N = (1−x 0 1 −x 0 4 )+(x 0 1 +x 0 4 ) exp(−δv 2 t 4s). The concurrence of this path exactly coincides with the upper concurrence bound derived in Eq. (8), which is not surprising because the distribution of concurrence shows sharp peaks along the concurrence bound. We note that the solutions for the low-concurrence branch (for v t = δv and v t = −δv) can be found numerically.
In the above analysis, for the half-parity case, a simple analytic solution for the most likely path for arbitrary values of v t was not forthcoming. However, if we further simplify the problem by considering a parity meter [6], then we can solve the equations of motion, Eqs. (10) and their conjugate equations, exactly. A parity meter has the same detector outputs for the even and odd parity subspaces, but the detector can distinguish between the subspaces. More detail of the calculation is presented in Appendix B.

C. Distribution of time to maximum concurrence
In the process of the entanglement generation, there are interesting quantities to investigate such as the maximum concurrence each individual trajectory can reach, and the time it takes to reach the highest value. We previously showed that qubit trajectories branch out to high and low concurrence subspaces. Therefore, one would expect that there are two most likely times for the qubit trajectories to reach their maximum concurrences (or their most entangled states).
We show in Figure 4 the normalized histograms of time for transmon qubit trajectories to reach their maximum concurrence. The histograms for the τ m = 0.36µs and 0.60µs measurement cases explicitly show double peaks, which agree with the branching of concurrence and the most likely qubit paths in Figure 3(b,c). The times at which these peaks are located can be theoretically predicted from the time-to-maximum-concurrence of the solutions of the most likely paths; as shown by the vertical dashed lines in Figure 3(b,c) and Figure 4: A 1,2 , B are the two most likely times to reach maximum concurrence (for low and high concurrence branches, respectively) for τ m = 0.60µs, and C 1,2 , D are the same but for the case with τ m = 0.36µs. The agreement between the theoretical prediction of the peaks and the peaks of the histograms are as good as the agreement of the theoryexperiment most likely paths in Figure 3(b,c). We note that for the weak measurement regime, τ m = 2.10µs, the bifurcation has not occurred yet during the measurement time T = 1.6µs. One would expect to see a branching effect, when the total measurement time is long enough.

V. CONCLUSION
We have investigated the process of entanglement generation between two spatially separated superconducting transmon qubits, and their statistical properties. The entanglement of the two qubits is created as a result of the half-parity dispersive measurement, via the microwave pulses sequentially interacting with both qubits. The strength of the joint measurement is arbitrary and we have studied three different values of the measurement strength. We examined the concurrence of individual trajectories and theoretically calculated its distribution from the quantum Bayesian approach, gradually projecting the two-qubit states to entangled states with high concurrence, and to separable states with zero concurrence.
The most likely path analysis was also carried out, predicting the most probable paths for the qubits trajectories. We found that in the two-qubit state space, there are three likely paths conforming to the three branches projecting to the 00⟩ subspace, the 11⟩ subspace, and the odd ( 01⟩, 10⟩) subspaces; the first two correspond to the lower branches of the concurrence bifurcation, and the last corresponds to the high concurrence branch. These theoretical most likely paths show excellent agreement with the experimental most likely ones extracted from the transmon trajectories data (for three independent data sets) based on the trace distance between trajectories in two-qubit state space. Moreover, we have presented the distributions of the time to the maximum concurrence for individual trajectories. The most likely path analysis was shown to be useful in predicting the peaks of these time distributions.
We conclude that the accurate tracking of quantum trajectories of a jointly measurement qubit system is possible, and that the physics of the entanglement creation statistics is well described by a quantum trajectory theoretical approach. The theoretical most likely paths to entanglement and concurrence distribution match the experiment excellently. This work shows the way to use this process as a control mechanism to entangle remote systems for quantum information processing purposes. In future work, similar questions can be posed about the full parity measurement, and we have made some predictions about that case in this work. Following the outline in Ref. [36,37], for the action principle for the continuous quantum measurement, the optimized paths starting from an initial state and ending at a final state after some time t is given by optimizing an action of the stochastic path integral, where the two-qubit variables x j and their conjugates p j for j = 1, 2, ..., 5 are implicitly functions of time t ′ . The first sum in the time integral is the logarithm of the joint probability density function of the measurement readout P({v t }) truncated to first order in dt, and the functional F j is the right hand side of Eq. (10). By extremizing the action Eq. (A1) over all variables x i , p i , v t , we get a set of 10 ordinary differential equations (ODEs) for the optimized path, and one equation for optimal measurement readout. The set of ODEs includes 5 differential equations of the two-qubit variables x i 's, as shown in Eq. (10), and 5 equations for the conjugate variables p i 's, x j A ij + x 5 p 5 B i + C i for i = 1, 2, 3, 4, (A2a) where, The optimal readout is given as a function of the twoqubit variables and the conjugate variables, connecting the ODEs for the qubit variables and the conjugate variables. By taking time-derivative of the function in Eq. (A3), and substituting both ∂ t x i and ∂ t p j with the ODEs above, we find that ∂ t v t = 0.