{'fraction_non_alphanumeric': 0.044196992308575364, 'fraction_numerical': 0.010255232847357748, 'mean_word_length': 3.9650422722523038, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 14, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Quantum Causal Unravelling

Complex processes often arise from sequences of simpler interactions involving a few particles at a time. These interactions, however, may not be directly accessible to experiments. Here we develop the first efficient method for unravelling the causal structure of the interactions in a multipartite quantum process, using only access to its input and output systems. Our method is based on a quantum algorithm whose complexity scales polynomially for all processes with bounded information loss. Under additional assumptions, we also provide a second algorithm that has lower complexity and requires only local state preparation and local measurements. Our algorithms can be used to identify processes that can be characterized efficiently with the technique of quantum process tomography. Similarly, they can be used to identify useful communication channels in a quantum networks, and to test the structure of uncharacterized quantum circuits.


I. INTRODUCTION
Many processes in nature arise from sequences of basic interactions, each involving a small number of physical systems. Determining the causal structure of these interactions is important both for basic science and for engineering. Often, however, the internal structure of a process of interest may not be directly accessible to experiments. For example, scattering experiments in high energy physics can probe the relation between a set of incoming particles and a set of outgoing particles, but typically cannot access the individual events taking place within the scattering region. In this and similar scenarios, a fundamental problem is to characterize the causal structure of the interactions by accessing only the inputs and outputs of the process of interest, while treating the intermediate steps as a black box. We call this problem, illustrated in Figure 1, the causal unravelling of an unknown physical process. Explicitly, the problem of causal unravelling is to determine whether an unknown process can be broken down into a sequence of simpler interactions, to determine the order of such interactions, and to determine which systems take part in each interaction. In principle, candidate answers can be extracted from a full tomographic characterization of the process under consideration. However, the complexity of process tomography grows exponentially in the number of inputs and outputs, making this approach unfeasible when the process involves a large number of systems.
In the classical domain, the problem of causal unravelling can be efficiently addressed with a variety of causal discovery algorithms [1][2][3]. Classical causal discovery algorithms often formulate causal relationships with graphical models and solve such models by structure learning algorithms, such as PC algorithm [1], Greedy Equivalence Search [4], and Max-Min Hill Climbing [5]. These algorithms cover a wide variety of problems, by making different sets of assumptions on the process under consideration. Typical assumptions include causal sufficiency-meaning that no variables are hidden-and causal faithfulness-meaning that the conditional independences among the variables are precisely those associated to an underlying graph used to model the causal structure. In general, however, causal discovery is intrinsically a hard problem: when no assumption is made, the complexity of all the known algorithms becomes exponential in the worst case over all possible instances [6,7].
In the quantum domain, the problem of causal unravelling is made even more challenging by the presence of corre- lations that elude a classical explanation [8,9]. In recent years, the quantum extension of the notion of causal model has been addressed in a series of works [10][11][12][13][14][15], providing a solid conceptual foundation to the field of quantum causal discovery. On the algorithmic side, however, the study of quantum causal models remained relatively underdeveloped. Specific instances of quantum causal discovery were studied in Refs. [16][17][18], showing that quantum resources offer appealing advantages. These examples, however, were limited to simple instances, typically involving a small number of variables and a small number of hypotheses on the causal structure. In more general scenarios, one approach could be to perform quantum process tomography and then to infer the causal structure from the full description of the process under consideration [19]. As in the classical case, however, the number of queries needed by a full process tomography grows exponentially with the number of systems involved in the process, making this approach impractical as the size of the problem increases.
In this paper, we provide an efficient algorithm for unravelling the causal structure of multipartite quantum processes without resorting to full process tomography. The algorithm is guaranteed to provide a correct answer in polynomial time whenever the process has bounded information loss, and satisfies a form of causal faithfulness property, namely that the strength of the causal relations, when present, is above a given threshold. When these assumptions are not satisfied, our algorithm can still be applied as a heuristic. Moreover, the efficiency can be further boosted in special cases, including i) the case where each input of a given interaction has a non-trivial causal influence on all the outputs of subsequent interactions, and ii) the case where the process belongs to special case of Markovian processes [12,19,20], where each output depends only on one previous input. In these cases, we devise an alternative algorithm that only requires local state preparations and local measurements, and the number of queries to the process is only logarithmic with the number of input and output wires.
The ability to infer the underlying causal structure of a process is useful for a variety of applications. Causal relations among quantum variables are relevant to the study of quantum networks [21][22][23], where the presence of a causal relation between two systems can be used to test whether it is possible to send signals from one node to another. In a realistic setting, the signalling patterns within a quantum network may be changing dynamically, depending on the number of users of the network at a given moment of time, on the way the messages are routed from the senders to the receivers, and also on changes in the environment, which may affect the availability of transmission paths between nodes. Such a dynamical structure occurs frequently in classical wireless networks [24,25], and is likely to arise in a future quantum internet. In this context, our algorithms provide an efficient way to detect dynamical changes in the availability of data transmission paths. Furthermore, characterizing the causal structure is useful even when the ultimate goal is to have a complete process tomography. Indeed, the tomography of a quantum process with a known causal structure can be significantly more efficient than the tomography of a generic process, owing to the fact that processes with a given causal structure are in one-to-one correspondence with a set of locally correlated multipartite states, for which there exist efficient tomographic protocols [26]. By first extracting information on the causal structure, it is therefore possible to choose the appropriate tomographic setup, sometimes achieving an exponential reduction of the number of samples.

II. RESULTS
Problem formulation. Let us start by giving a precise formulation of the problem of quantum causal unravelling. In this problem, an experimenter is given access to a multipartite quantum process, with inputs labelled as A 1 , . . . , A nin , and outputs labelled as B 1 , . . . , B nout . Note that, in general, different labels may refer to the same physical system: for example, system A 1 could be a single photon with a given frequency, entering in the interaction region, and system B 1 could be a single photon with the same frequency, exiting the interaction region. Mathematically, the process is described by a quantum channel, that is, a completely positive trace-preserving (CPTP) linear map C transforming operators on the tensor product space H A1 ⊗ · · · ⊗ H An in to operators on the tensor product space H B1 ⊗ · · · ⊗ H Bn out . Note that, without loss of generality, one can always assume n in = n out = n, as this condition can be satisfied by adding a number of trivial systems with one-dimensional Hilbert space. In the following, we will denote by L(H) be the set of linear operators on a generic Hilbert space H, and by S(H) be the subset of density operators on H.
Decomposition of a physical process into a sequence of k interactions. The inputs (outputs) of the processes are divided into k non-overlapping subsets. The i-th subset of the inputs (outputs) consists of the systems that enter (exit) the interaction at the i-th step (rectangular boxes in the picture). The wires connecting one interaction to the next represent intermediate systems that in this stage are not directly accessible to experiments. For example, a photon could enter the first interaction, remain as an intermediate system between the first and second interaction, and then exit the interaction region after the second interaction.
The problem of causal unravelling is to determine whether a multipartite process can be broken down into a sequence of interactions, as in Figure 2, and, in the affirmative case, to determine which systems are involved in each interaction. Mathematically, the problem is to find a partition of the set {B 1 , . . . , B n }, such that the multipartite process can be decomposed into a sequence of interactions, with the i-th interaction involving input systems in P i and output systems in Q i . A compact description of such sequences is provided by the framework of quantum combs [27,28], which provides a set of necessary and sufficient conditions for determining whether a given process admits a causal unravelling with partitions {P i } and {Q i }. We say a process C has a causal unravelling (P 1 , Q 1 ), . . . , (P k , Q k ) if it can be decomposed into the form of a quantum comb with k teeth as in Figure 2. . This is a bipartite graph because inputs can be independently controlled, and every two outputs are independent conditioned on all inputs.
The causal unravelling of a given multipartite process determines the possible signalling relations between inputs and outputs. With respect to this decomposition, input systems at a given time can only signal to output systems at later times. The resulting pattern can be graphically illustrated by a causal graph [1,2,14], as shown in Figure 3. Note that the causal graph includes all the signalling relations that are in principle compatible with the structure of the interactions. However, a quantum specific process may not exhibit any signalling from a specific input to a specific output, even though the causal structure of the interactions would in principle permit it. When this occurs, the causal unravelling of a process may not be unique: for example, a bipartite process with inputs {A 1 , A 2 } and outputs {B 1 , B 2 }, with the property that A 1 cannot signal to B 2 and A 2 cannot signal to B 1 admits two causal unravellings: one where A 1 and B 1 appear before A 2 and B 2 , and one where A 2 and B 2 appear before A 1 and B 1 . In this case, any of the two options is a valid solution of the causal unravelling problem.
In its general form, the causal unravelling problem is hard, even in the classical setting. However, we will now show that, under a few assumptions, the quantum version of the problem can be solved efficiently. The first assumption is that the basic interactions appearing in the causal unravelling involve a small number of systems, independent of the number of inputs and outputs. At the fundamental level, this assumption is motivated by the fact that interactions are local, and typically involve a small number of systems. Mathematically, the assumption is that the cardinality of all the sets in the partitions is no larger than a constant c independent of n. For simplicity, we will first restrict our attention to the special case c = 1, meaning that the process can be broken down into interactions involving only one input and one output at a time. Hereafter, we will denote by Comb[(A 1 , B 1 ), . . . , (A n , B n )] the set of quantum combs with n teeth where the i-th tooth has input A i and output B i .
To quantify the efficiency of our algorithms, we will focus on the sample complexity, namely the number of blackbox queries to the channel C, and on the computational complexity, including additional quantum and classical computation time measured by the number of elementary quantum gates and classical operations.
Efficient quantum causal unravelling. We now provide an efficient quantum algorithm for quantum causal unravelling. The main idea of the algorithm is to recursively find the last interaction in the decomposition of a given process. To illustrate this idea, consider the example of Figure 3. When we remove B 4 , the node A 4 becomes disconnected from all the other nodes, meaning that the state of system A 4 does not affect the state of the other systems. By testing this independence condition, we can in principle check whether the graph of Figure 3 is an appropriate model for the process under consideration. In general, suppose (A x , B y ) are the input and output of the last interaction. Since A x signals to only B y , if we ignore B y , A x is independent with the joint system formed by all systems excluding A x and B y . This gives the following criterion that the last interaction must satisfy: Proposition 1. [27] The last interaction of a process C involves the input/output pair (A x , B y ) if and only if A x is independent from A =x B =y , the joint system containing all systems other than A x and B y .
More details can be found in Section A of the Appendix. By scanning the possible pairs (A x , B y ), we can find if one of them satisfy the above criterion, and in the affirmative case, we can assign that pair to the last interaction. Note that, in general, there may be more than one pair that satisfy the required condition. After one pair is found, one can reduce the problem to a smaller graph containing n − 1 inputs and n − 1 outputs. If at every step a suitable pair is found, then the final result is a valid causal unravelling of the original process. If at one step no pair can be found, the algorithm will then conclude that no causal unravelling with c = 1 exists for the given last pair. At this point, the algorithm can continue by considering causal unravellings with higher values of c for the remaining subgraph.
A description of the algorithm is provided in Algorithm 1. The algorithm runs in a recursive manner: for an n-tooth comb, it finds the last tooth of the comb, remove the tooth (feeds the input with an arbitrary state and discards its output), and reduces the problem to finding the causal unravelling of an (n − 1)-tooth comb. We repeat the above procedure until we reach the bottom case n = 1, and thus obtain the order of all inputs and outputs. If the last tooth cannot be found in some iteration, it means that the current channel cannot be further decomposed, and the algorithm will output the trivial causal unravelling (P, Q) where P (Q) is the set of all input (output) wires of the current channel. Algorithm 1: Efficient quantum causal unravelling algorithm Input : Black-box access to quantum channel C Output : A causal unravelling of C 1 if C has only one input wire A x and one output wire B y then 2 Output (A x , B y ) and exit this function; is the last tooth then // Done by checking whether A x and A =x B =y are independent using the quantum circuit described in the Methods 6 Let C A =x ,B =y (ρ) := Tr By [C(ρ ⊗ τ Ax )], where τ Ax is an arbitrary state on system A x ; // Remove A x and B y from C 7 Recursively run this algorithm on the reduced channel C A =x ,B =y ; Let T test (n) be the number of independence tests required for an n-to-n channel. In this algorithm, an independence test is performed for at most each of the n 2 input-output pairs (A x , B y ), resulting in O(n 2 ) independence tests. After finding a last tooth, the problem size is reduced to n − 1. Therefore, T test (n) can be given by the recursive relation To show the efficiency of this algorithm, the remaining step is to show that the independence tests can be efficiently realized. This is a non-trivial problem, because testing whether two generic systems are in a product state is computationally hard in the worst-case scenario [29]. Nevertheless, in the Methods we design a quantum circuit that performs the independence tests efficiently under assumptions of bounded information loss and causal faithfulness. Our circuit converts the independence test to the estimation of the distance between quantum states, which is done by the SWAP test [30]. Now, we give the efficiency guarantee for our algorithm. An important parameter entering into the analysis is the degree of independence between systems, defined in the following. For a state ρ, we say two disjoint subsystems S and T are independent if ρ ST = ρ S ⊗ ρ T , where ρ S , ρ T and ρ ST are the marginal states of ρ on the subsystems S, T and the joint system ST , respectively. The degree of independence between S and T is then defined as χ 1 (S; T ) ρ := X † X denotes the trace norm. Clearly, χ 1 (S; T ) ρ = 0 if and only if ρ ST = ρ S ⊗ ρ T . More generally, the trace distance is related to the probability to distinguish the states, and thus χ 1 (S; T ) ρ measures the probability that an observer correctly decides whether the subsystems are independent or correlated. In the following, we will apply this definition to the Choi state of the process C, defined as the state C := (C ⊗ I)(|I I|)/d in , where d in is the total dimension of all the input systems, and |I = din j=1 |j ⊗ |j is the canonical (unnormalized) maximally entangled state. With this choice, χ 1 (S; T ) C satisfies the conditions for a quantum causality measure, as defined in Ref. [31].
To facilitate the efficiency analysis of our algorithm, we first put the process C in a standard form where all input and output wires have the same dimension d A . This standard form does not limit the generality of the quantum process we investigate. For a process C with input dimensions d A1 , . . . , d An and output dimensions d B1 , . . . , d Bn , we can pick d A = max{d A1 , . . . , d An , d B1 , . . . , d Bn } and regard each input or output wire of C as a subspace of d A -dimensional system, thus transforming C into a process whose wires all have dimension d A . In our analysis, we will use the Kraus rank of the process in the standard form to characterize the information loss. We assume that the information loss is bounded, which is given by the following assumption: Assumption 1. The Kraus rank of C, after transforming it into the standard form, is bounded by a polynomial of n.
To ensure that the algorithm outputs the correct causal unravelling, we further require that the process satisfies a form of causal faithfulness, meaning that the strength of causal relations are either zero or above a threshold. Using χ 1 as a quantitative measure of correlation, we adopt the following assumption: There exits a number χ min > 0 such that, for any two disjoint sets of wires S and T being tested for independence, either where C is the Choi state of the quantum process C.
The threshold χ min determines the resolution of the independence tests. To guarantee the correctness of Algorithm 1, the independence tests must be precise enough to detect correlations above this threshold with high probability. The efficiency and correctness of Algorithm 1 are given in the following theorem, whose proof is in the Appendix: Theorem 1. Under Assumptions 1 and 2, for any confidence parameter κ 0 > 0, Algorithm 1 satisfies the following conditions: 1. With probability 1 − κ 0 , the output of Algorithm 1 is a correct causal unravelling for C.
2. The number of queries to C is in the order of where r C is the Kraus rank of C in the standard form with all wires having dimension d A .
3. The computational complexity is in the order of O(T sample n log d A ).
Theorem 1 guarantees that, under appropriate assumptions, the sample complexity of our algorithm is polynomial in the number of input and output systems of the process under consideration. This feature is in stark contrast with the exponential complexity of full process tomography. As a consequence, our algorithm offers a speedup over for other algorithms, such as the one proposed in Ref. [19], which require process tomography as an intermediate step.
On the other hand, if no assumptions on the process are made, it is known that the sample complexity is exponential, even in the classical case.
Causal unravelling with local observations. We now show that, under some assumptions on the input-output relations, one can design algorithms that have much lower sample complexity compared to Algorithm 1, and are more experimentally friendly, in that they require only local state preparation and local measurements.
In the Methods, we show an efficient algorithm to detect the pairwise correlations between input and output wires of C with local state preparation and local measurements. The algorithm computes a Boolean matrix ind ij such that, with high probability, for every i and j, A i and B j are approximately independent whenever ind ij = true, and are correlated whenever ind ij = false. With some assumptions, this Boolean matrix ind ij is sufficient to give the exact causal unravelling. The first case is given by the following assumption on the process C: (1) ), . . . , (A σ(n) , B π(n) )], and there exists a constant χ min > 0 such that, for any pair of input and output wires A σ(i) and B π(j) , if j ≥ i, then χ 1 A σ(i) ; B π(j) C ≥ χ min .
Assumption 3 indicates a non-trivial correlation between any input-output pair that the input is before the output. In other words, for any j ≥ i, C A σ(i) ,B π(j) is away from C A σ(i) ⊗ C B π(j) by distance χ min . Meanwhile, Assumption 3 defines a total order of the input (output) wires, and ensures a unique causal unravelling that C is compatible with. Under this assumption, if A i is the k-th input, namely i = σ(k), A i is correlated with n − k + 1 output wires including every output B π(j) with j ≥ k, and is independent with the other outputs. In other words, if we find A i is correlated with exactly c A (i) output wires, it must be the (n − c A (i) + 1)-th input. With this, the order of input wires can be exactly determined, and a similar statement can be applied to order the output wires.
In Section B of the Appendix, we give the details of this algorithm, and analyze its efficiency given by the following theorem: satisfying Assumption 3, there is an algorithm that satisfies the following conditions: 1. With probability 1 − κ, the algorithm outputs the correct causal unravelling (A σ(1) , B π(1) ), . . . , (A σ(n) , B π(n) ).
2. The algorithm uses only local state preparations and local measurements and the number of queries to C is in the order of where 3. The computational complexity is in the order of O(N n(n + d A + d 4 B )).
Note the sample complexity of this algorithm grows only logarithmically with n.
A similar idea could be adopted to the case where the process belongs to a special case of Markovian processes [12,19,20], where each output depends only on one previous input. This indicates that the process is decomposable to a tensor product of n channels each with one input and one output. In our problem, the order of inputs and outputs is unknown, and we have the following assumption: Since each output is related to at most one input and each input affects at most one output, after obtaining ind ij , we can obtain the causal unravelling by matching each input-output pair (A i , B j ) with ind ij = false.
If there exists a threshold χ min > 0 such that either χ 1 (A i ; B j ) C = 0 or χ 1 (A i ; B j ) C ≥ χ min holds for every A i and B j , then the algorithm has the same complexity as in Theorem 2 that is logarithmic in n. However, in this specific case, we give an additional efficiency guarantee with an error bound defined by the diamond norm in the following theorem, whose proof is in the Appendix.
Theorem 3. For a quantum process C satisfying Assumption 4, there is an algorithm that outputs a causal unravelling (A 1 , B π(1) ), . . . , (A n , B π(n) ) satisfying the following conditions: 1. With probability 1 − κ, the causal unravelling is approximately correct in the following sense: 2. The algorithm uses only local state preparations and local measurements and the number of queries to C is in the order of 3. The computational complexity is in the order of O(N n(n + d A + d 4 B )). We have used the the diamond norm [32], also known as completely bounded trace norm, as a distance measure between channels, which is defined for C, D : Causal unravelling with interactions between more inputs and outputs. In the algorithms shown so far, we assumed that the process under consideration admits a causal unravelling where each interaction involves exactly one input and one output. More generally, Algorithm 1 can be easily extended to the scenario where each interaction involves at most c inputs and c outputs of the original process. Instead of considering each wire separately, the idea is to consider a subset of at most c input (output) wires and perform independence tests on the subsets.
In Algorithm 1, one enumerates an input-output pair (A x , B y ) and checks whether it is the last tooth by performing an independence test between A x and A =x B =y . To deal with larger c, we replace this procedure by enumerating a subset of input wires P ⊂ {A 1 , . . . , A nin } and a subset of output wires Q ⊂ {B 1 , . . . , B nout }, satisfying |P | ≤ c and |Q| ≤ c. Then we check whether (P, Q) is the last tooth of C, which, according to Proposition 1, is equivalent to checking the independence between P and A / ∈P B / ∈Q : The independence tests can still be implemented with SWAP tests. Like Algorithm 1, after we decide (P, Q) to be the last tooth, the wires in P and Q are removed from consideration, and the problem is reduced to the causal unravelling of a smaller channel. This process is done recursively until one reaches the bottom case. We give the detailed algorithm and analysis in Section F of the Appendix. For constant c, under some assumptions, the complexity of the algorithm is still polynomial in n and d A , while the exponent depends on c.

III. DISCUSSION
In this work we provided an algorithm for quantum causal unravelling, with polynomial sample complexity whenever the process under consideration has a bounded Kraus rank. For special cases, we also proposed more efficient algorithms that require only local state preparation and local measurements. All of these algorithms are guaranteed to be correct under some assumptions on the process, including the assumption that the causal relations between variables are above a certain threshold. When these assumptions are not satisfied, the algorithms can still be used as a heuristic, and may still produce approximately correct solutions, in the sense that they may just fail to detect causal relations that are too weak to be practically relevant. On the other hand, the lack of a universally valid algorithm for identifying all causal relations is a well known issue also in the classical case, where most existing algorithms take exponential time in the worst case.
Our solution to the quantum causal unravelling problem discovers a linear structure of the process given by a sequence of input and output wires. This is a partial solution to the more general quantum causal discovery problem, whose goal is to produce a full causal graph describing all the causal correlations between systems. Our algorithm can be used as the first step for quantum causal discovery, and to obtain the full causal structure, additional tests maybe adopted to detect signalling between more subsets of input and output wires. Since the most general quantum causal discovery problem is intrinsically hard, it will be demanding to find out to what extent, beyond the linear structure produced in this work, the problem can be solved by an efficient algorithm.
Besides the identification of causal relations per se, one of the applications of our algorithms is to identify accessible communication channels in quantum networks. The assumption of a minimal strength of the causal relations fits well this application, since the strength of a causal relation between two systems is related with the amount of information that can be transmitted from one system to the other. Note that the mere dependence of a quantum system from another enables the transmission of a non-zero amount of classical information, and is a necessary condition for the transmission of quantum information. In this respect, our algorithms can serve as a first test of dependency, which can be used to decide whether a full process tomography is worth being performed.
A related application is to process tomography itself. Our causal unravelling algorithms can be used to detect the causal structure of multipartite quantum processes, and knowing the causal structure of the system may allow us to ignore unnecessary correlations. Thanks to this pre-tomography step, one can find out that the given process can be efficiently represented, for example, by a tensor network [33,34], for which tomography can be performed efficiently [26]. Based on the causal structure and tomography data, we may model the underlying physical process by finding the smallest quantum model that reproduces the observed data [35][36][37]. More generally, characterizing the causal structure of a multipartite process as a tensor network enables the use of efficient protocols that exploits the tensor network structure, including simulation protocols [34,38,39] and compression protocols [40].

IV. METHODS
Efficient tests for the last tooth via the SWAP test. In the Results, we have given the framework of the algorithm. In this section, we discuss how the tests for the last tooth, namely line 5 of Algorithm 1, can be carried out efficiently.
Consider a process C of three input wires and three output wires as in Figure 4, and suppose that we want to test whether (A 1 , B 1 ) is the last tooth, which is equivalent to the independence test between A 1 and A 2 A 3 B 2 B 3 . To test the independence between two sets of wires of a quantum process, it is convenient if one set contains only input wires and the other contains only output wires. By adding maximally entangled states, we can make the conversion as shown in Figure 4.
denotes maximally entangled states on A2A 2 and A3A 3 . By this conversion, the wire in the dotted frame and wires in the dashed frame become the input and outputs of the new channel CA 1 →A 2 A 3 B 2 B 3 .
Testing the independence between A 1 and A 2 A 3 B 2 B 3 is equivalent to checking whether the quantum channel C A1→A2A3B2B3 defined in Figure 4 is a constant channel or not. To test this, we simply observe whether the outputs of C A1→A2A3B2B3 depend on the inputs.
To check whether {ρ α } are all the same, we pick ρ 1 as the pivot, and estimate its distances to ρ α with α > 1. We assert C A1→A2A3B2B3 to be a constant channel if all estimated distances are below some threshold δ.
The distances between ρ 1 and ρ α can be estimated via SWAP tests [30]. The SWAP test uses the quantum circuit in Figure 5 to estimate Tr[ρσ] for two given quantum states ρ and σ. If we run the circuit in Figure 5 for N times and let c + be the number of times observing outcome |+ , then 2c + /N − 1 is an estimate of Tr[ρσ]. The algorithm that yields an estimate of Tr[ρσ] is as follows, which produces an estimate with error no more than ε with probability 1 − κ as shown in Lemma 1. Proof. For each run, the probability that the outcome is |+ is (1 + Tr[ρσ])/2. By Hoeffding's inequality, The SWAP test is efficient in the sense that its circuit complexity is linear in the number of qubits representing the systems. For d-dimensional states ρ and σ, the controlled-SWAP gate acting on ρ and σ can be realized with log d controlled-SWAP gates acting on each pair of corresponding qubits of ρ and σ, and thus can be implemented with O(log d) gates.
To summarize, our method uses the equivalence between the following statements: 1. (A x , B y ) is the last tooth of C. 3. The channel C Ax→A =x B =y as shown in Figure 6 is a constant channel.

A x and
4. Pick an informationally complete set of states {ψ α } for wire A x . Let ρ α := C Ax→A =x B =y (ψ α ) be the output of C Ax→A =x B =y . Then for any α, ρ 1 − ρ α 2 is below some threshold.
Based on the SWAP test, the algorithm to test for the last tooth is as follows. The choices of the threshold δ and confidence κ are given during the proof of Theorem 1 in Section B 2 of the Appendix. Given the subroutine to check the last tooth, the main recursive algorithm is shown in Algorithm 1.
Testing the independence between all input-output pairs. In this section, we show the algorithm to test the independence between every pair of input and output wires, resulting in a Boolean matrix ind ij . The first step is to collect statistics of the channel C. When collecting data, we need to ensure that the information encoded in the quantum data are preserved, for which purpose we use informationally complete POVMs [41]. The elements of such a POVM on Hilbert space H form a spanning set of L(H), and the number of elements can be chosen as (dim H) 2 .
We pick an informationally complete POVM {P Ai,α }  β=1 for each output wire B j . We pick N to be the number of queries to the channel C. In each query, we measure the Choi state of C with the POVM {P A1,α1 ⊗ · · · ⊗ P An,αn ⊗ Q B1,β1 ⊗ · · · ⊗ Q Bn,βn } α1,...,αn,β1,...,βn This POVM can be realized with local measurements on each wire. Let (α n ) be the outcome of k-th query. We can write the outcomes in a matrix as shown in Figure 7.
FIG. 7. Matrix of measurement outcomes. Each row corresponds to a query to the process, and each column corresponds to an input or output wire. When we investigate Ai and Bj, we only look at the corresponding columns in the matrix and check the independence from data in these columns.
This matrix will be used to determine the independence between every pair of input and output wires. When we investigate A i and B j , we only look at the corresponding columns in the matrix and check the independence from data in these columns. Since we are using an informationally complete POVM, the independence between the two columns of the matrix is equivalent to the independence between the input A i and output B j . With this idea, we will be able to compute a Boolean matrix ind ij such that ind ij = true if and only if A i and B j are independent, up to an error related to N , the number of queries.
One may notice that the POVM in Eq. (7) is enough for a full process tomography of C. This is true if we pick an exponentially large N as large as the square of the product of all input and output dimensions. However, to compute ind ij with a small error, N does not need to be large. We show that N could be chosen to grow only logarithmically with n, far less than a full process tomography.

Lemma 2.
Given a channel C with input wires A 1 , . . . , A n and output wires B 1 , . . . , B n , there exists an algorithm that satisfies the following conditions: 1. With probability 1 − n 2 κ 0 , the algorithm produces an output satisfying: where χ − > 0 is a threshold that can be freely chosen.

The number of queries to C is in the order of
The detailed algorithm and the proof are in Section C of the Appendix. The Boolean matrix ind ij gives a partial order of the wires: if ind ij = false, then A i must be before B j . This partial order may not produce the full ordering of the wires, but will be a convenient initial guess for the causal unravelling, since its complexity is much lower than the general algorithm Algorithm 1. In the Results, we show that under Assumptions 3 or 4, this partial order is enough to infer the full order.
Furthermore, the matrix of outcomes (Figure 7) is a conversion from the quantum process to classical data, making it possible to adopt classical causal discovery algorithms [1,[3][4][5]. For example, one could use algorithms based on the graph model, and try to find a graph in the form of Figure 3 that best fits the observation data in Figure 7.
In this section, we give the theoretical background of our algorithms. In the main text, we showed that independence tests can be used to find the last tooth of a comb. Here we justify this idea by the following proposition that gives a sufficient and necessary condition for (A x , B y ) to be the last tooth. The condition is based on the Choi state [42] of a channel, defined as C ∈ H A1 ⊗ · · · ⊗ H An ⊗ H B1 ⊗ · · · ⊗ H Bn , ..An is a maximally entangled state, is a basis of H A1 ⊗ · · · ⊗ H An . The condition Eq. (A1) for the Choi state C is equivalent to that for the process C, the sets of systems {A x } and {A i |i = x} ∪ {B j |j = y} are independent.
When we collect classical statistics from a quantum state or process, we need to ensure that the information encoded in the quantum data are preserved, for which purpose we introduce the informationally complete POVM [41].

Definition 1.
A POVM {P α } on Hilbert space H is informationally complete if its elements form a spanning set of linear operators on H, namely for any X ∈ L(H), there exist complex numbers {p α } such that X = α p α P α .
For efficiency considerations, we prefer to choose an informationally complete POVM with the minimal number of elements. To form a spanning set of L(H), the minimal informationally complete POVM contains dim(H) 2 elements, which are linearly independent. A typical example of such POVMs is the symmetric informationally complete POVM (SIC-POVM) [43], which has been found for most relevant low-dimensional Hilbert spaces.
The frame operators defined below will be used to characterize the sensitivity of measurement outcomes about the state being measured.
A frame operator with a larger minimum eigenvalue indicates that the outcome probabilities are more sensitive to any change of the state being measured, and the POVM is more efficient for obtaining the full information of the state. We will see the consequences of the minimum eigenvalue of frame operators in the analysis of our algorithms.
In the following sections, we will use the following inequality between the trace distance and the Hilbert-Schmidt distance [44]: Appendix B: Proof of Theorem 1 We prove Theorem 1 in this section. We assume that the process C is already in its standard form with all input and output wires having dimension d A , and thus r C = rank(C).
Independence test. The coefficients {p x } in Definition 1 can be computed with the following lemma: For an informationally complete POVM {P x }, let F := x |P x P x | be its frame operator. Then any linear operator X ∈ L(H) can be decomposed as X = x p x P x where p x = P x | F −1 |X . The same is true for an informationally complete set of states {|ψ x }: Proof.
where {|P x } forms a basis of space H ⊗ H, and p x are the coefficients of |X in this basis. From we obtain the expression of p x . The expression of p x is identical by substituting P x with ψ x .
Since {P x } forms a basis for linear operators in H, {|P x } forms a basis for H ⊗ H, and F is always full-rank. Recall that to use the SWAP test for independence test, we need the equivalence between the following statements: 1. A x and A =x B =y are independent.
2. The channel C Ax→A =x B =y as shown in Figure 6 is a constant channel.
3. Pick an informationally complete set of states {ψ α } for wire A x . Let ρ α := C Ax→A =x B =y (ψ α ) be the output of C Ax→A =x B =y . Then for any α, ρ 1 − ρ α 2 is below some threshold.
where D is the Choi state of D, and λ min is the minimum eigenvalue of the frame operator of Proof. The diamond norm can be evaluated as: where |Φ is a pure state on systems A and A , and A is a reference system with dim is an informationally complete set, by Lemma 3, one has where p ijk = ψ k | F −1 |α i |α j with F := k |ψ k ψ k | being the frame operator. Then Note that ∀i, j, Last, be an informationally complete set of states for system A. Let ρ α := C(ψ α ). If ∀α, ρ α − ρ 1 1 ≤ ε, then where λ min is the minimum eigenvalue of the frame operator of {ψ k } Then be an informationally complete set of states for system A. Let ρ α := C(ψ α ). If ∀α, ρ α − ρ 1 2 ≤ ε, then where λ min is the minimum eigenvalue of the frame operator of {ψ α } d 2 A α=1 , and rank(C) is the Kraus rank of C. Proof. According to Eq. (A2), we have ρ α − ρ 1 2 1 ≤ (rank(ρ α ) + rank(ρ 1 )) ρ α − ρ 1 2 2 ≤ 2rank(C)ε 2 = 2rank(C)ε 2 . Substituting this into Lemma 5 we get Eq. (B39).

Bounding the rank during recursion
In the following, we first give some lemmas on the Kraus rank of channels, and then prove Lemma 8, which can be used in the inductive proof of Theorem 1 to upper bound the rank of C in any recursion.
Lemma 7. For a constant channel C : S(H A ) → S(H B ) satisfying C(X) = Tr[X]σ B , one has rank(σ B ) = rank(C)/ dim H A , where rank(C) is the Kraus rank of channel C.
Proof. Since C(X) = Tr[X]σ B , the Choi state of C is C = I A /d A ⊗ σ B . Then rank(C) = rank(I A )rank(σ B ), and rank(σ B ) = rank(C)/rank(I A ) = rank(C)/ dim H A . Then Lemma 8. Let C (n) be the initial value of C in Algorithm 1, and let C (k) , which has k input wires and k output wires, be the value of C in the (n − k + 1)-th recursive call of Algorithm 1. Define the follow propositions: • P (k): in the (n − k + 1)-th recursive call of Algorithm 1, the pair (A x , B y ) passing the test is indeed the last tooth, namely • Q(k): rank(C (k) ) ≤ rank(C (n) ).
Proof. According to P (k), C (k) ∈ Comb[(A =x , B =y ), (A x , B y )], and the following channel C (k) Ax→A =x B =y is shown in Figure 6.

Proof of Theorem 1
Now we are ready to prove Theorem 1.
Proof of Theorem 1. The SWAP test produces an estimate within error ε with probability no less than 1−κ according to Lemma 1. There are no more than 2n 3 d 2 A SWAP tests, and the probability that all of them are within error ε is no less than 1 − 2n 3 d 2 A κ. Let κ = κ 0 /2n 3 d 2 A , and from now we assume all SWAP tests are within error ε, which has probability no less than 1 − κ 0 .
Let C (n) be the initial value of C in Algorithm 1, and let C (k) , which has k input wires and k output wires, be the value of C in the (n − k + 1)-th recursive call of Algorithm 1.
Consider testing whether (A x , B y ) is the last comb of C (k) with Function 2. If (A x , B y ) passes the test, we know that for all α, p 1 + p α − 2p 1α ≤ δ. Since p 1 + p α − 2p 1α is an estimate for Tr[ρ up to error 4ε, we infer that ρ 1 − ρ α 2 2 ≤ δ + 4ε.
Let ε = δ/4, we have ρ 1 − ρ α 2 ≤ √ 2δ. Let S be the remaining input wires excluding A x , and T be the remaining output wires excluding B y . By Lemma 6, min rank(C (k) )δ. If we set the threshold δ ≤ χ 2 min λmin 17d 5 A rank(C (k) ) , then we have χ 1 (A x ; S ∪ T ) C (k) < χ min . Since χ 1 (A x ; S ∪ T ) C (k) = χ 1 (A x ; S ∪ T ) C (n) , by Assumption 2, χ 1 (A x ; S ∪ T ) C (k) = 0. Therefore, if (A x , B y ) passes the test of Function 2, then (A x , B y ) must be the last tooth of the comb, in some causal unravelling of C (k) . Using the notation in Lemma 8, P (k) holds for this iteration.
Take the definition of P (k) and Q(k) from Lemma 8. Let r C := rank(C (n) ). If Q(k) holds, we have rank(C (k) ) ≤ r C . In any recursive call of Algorithm 1, if Q(k) holds, then we can choose δ = A rank(C (k) ) , and by the argument above, P (k) holds. Combining this with Lemma 8, we have Q(k) ⇒ P (k) and Q(k) ∧ P (k) ⇒ Q(k − 1), meaning Q(k) ⇒ Q(k − 1). Last, since Q(n) obviously holds according to Assumption 1, we conclude that Q(k) and P (k) holds for every k, showing the correctness of our algorithm.
According to Function 1, The sample complexity is then and the computational complexity is O(T sample n log d A ) from the realization of SWAP tests. Last, recall that λ min is the minimum eigenvalue of the frame operator of an informationally complete set of states

Appendix C: Independence test algorithm between all input-output pairs
Our independence test algorithm will be based on the estimation of χ 1 (ρ A,B ). The idea is essentially quantum state tomography: measure ρ AB with informationally complete POVMs {P α } on system A and {Q β } on system B, use the outcome statistics to reconstruct ρ AB , and estimate χ 1 (ρ A,B ) based on the reconstructed state. The algorithm is given in Function 3.
Lemma 9. For any positive real number ε > 0, over the N measurement outcomes on independent copies of ρ AB , with probability 1 − κ 0 , the output of Function 3 estimates χ 1 (ρ A,B ) with the following error bound That is, the inequality holds for any positive real number ε > 0. Here, Pr(C) expresses the probability that the condition C is satisfied.
In other words, to reach confidence 1 − κ 0 and error ε 0 , the number of samples N is of order [43]. In this case, the number of samples Function 3 could be used on quantum channels to discover the correlation between input and output wires. To characterize a quantum channel, we could transform it to its Choi state by feeding the input with half of a maximally entangled state, and apply Function 3 on the Choi state. To save the cost of creating and maintaining entangled states, we could use Function 3 in an alternative way. Let {ψ α } β=1 . This process generates the same statistics for (α, β) according to the following equation: is the Choi state of channel C. The left hand side of Eq. (C5) is the probability of sampling ψ α times the probability of outcome β on the output conditioned on the input being ψ α , and the right hand side of Eq. (C5) is the probability of outcome (α, β) if one measure the Choi state C with POVM {P α ⊗ Q β }. Now we adopt Function 3 to causal unravelling. By testing the correlation between every pair of input and output wires, one is able to infer the causal unravellings: if output B j is correlated to input A i , B j must be after A i . In this section, we study cases where such reasoning yields the full causal unravelling.
We introduce an algorithm that discovers the correlations between every pair of input and output wires. The idea is to first collect all the required classical data, and deduce all the correlations according to the statistics. During this procedure, the data are reused for different input-output pairs, leading to fewer queries to the channel C.  To understand Function 4, we observe that {(a are sampled from the same distribution as the outcomes of measuring C Ai,Bj with {P i α ⊗ Q j β }, and estimatechi1 produces an estimate of χ 1 (A i ; B j ) C . Conditioned on all n 2 number of calls to estimatechi1 succeed, which has probability no less than 1 − n 2 κ 0 , all the estimations χ 1 (A i ; B j ) C have error no greater than ε 0 . We summarize these observations in the following lemma: Lemma 10. With probability 1 − n 2 κ 0 , Function 4 produces an output satisfying: where κ 0 and ε 0 are related as , λ min,in := min λ min (F i ), λ min,in := min λ min (G i ), and F i (G i ) is the frame opertor Lemma 11. Let {P α } be an informationally complete POVM and F := α |P α P α | be its frame operator. For two Hermitian operators ρ and σ, define p α := Tr[P α ρ] and q α := Tr[P α σ]. Then one has where X 2 := Tr[X † X] denotes the Schatten 2-norm, also known as the Frobenius norm, and λ max (F ) and λ min (F ) are the maximum and minimum eigenvalues of F , respectively.
The running time of Algorithm 2 mostly attributes to Function 4. We first discuss how to choose N , the number of queries to C. The relationship between N , ε 0 and κ 0 is given in Lemma . By our choice of parameters ε 0 = χ min /3 and κ 0 = κ/n 2 , we have To analyze the computational complexity, we consider the two main parts of Function 4: the query part and the estimation part. The query part invokes the channel N times, and each time n state preparation and n measurements are performed. Here we assume that state preparation on a d-dimensional Hilbert space takes time O(d) [45] and performing a POVM with m outcomes takes time O(m 2 ), according to the Naimark's dilation theorem that converts m-outcome POVM to projective measurement via a m-dimensional unitary gate [46] implementable with O(m 2 ) basic gates [47]. To perform the POVM with m = where C Ai→B π(i) (ρ Ai ) := Tr B =π(i) [C(ρ Ai ⊗ I A =i /d A =i )] is channel CC restricted on input A i and output B π(i) , and C Ai B π(i) (ρ Ai ) := Tr[ρ Ai ]C B π(i) is a constant channel that outputs C B π(i) , the marginal Choi state of C on system B π(i) . Clearly, D ∈ Comb[(A 1 , B π(1) ), . . . , (A n , B π(n) )], since it is a tensor product of channels from A i to B π(i) . Now we show that the difference between C and D is small.
According to Assumption 4, C = n i=1 C Ai→B π (i) = i∈A matched C Ai→B π (i) ⊗ i / ∈A matched C Ai→B π (i) . For channel C and each i ∈ A matched , A i is correlated to B π(i) but not correlated to all other outputs, and by the way π is computed, one must have π(i) = π (i). Therefore, and the difference between C and D only occurs in the second part involving i / ∈ A matched . Since is the constant channel defined as C Ai B π (i) (ρ Ai ) := Tr[ρ Ai ]C B π (i) . Eq. (E3) then becomes Consider the channel C Ai→B π (i) . For i / ∈ A matched , ind i,π (i) = true and with high probability χ 1 A i ; B π (i) C ≤ 2χ − . By definition of χ 1 , Eq. (E6) is the distance between the Choi states of C Ai→B π (i) and C Ai B π (i) , which gives a bound of the diamondnorm distance between the channels: Combining Eqs. (E5) and (E7), we obtain where d A := max i d Ai . The sample complexity of Algorithm 3 is given by Lemma 2 as Let κ := n 2 κ 0 and ε := 2nd A ε 0 , we have and with probability 1 − κ, C − D ≤ ε. The running time of Algorithm 3 is the same as Algorithm 2 and has been discussed in Section D.
complete set {ψ α } d 2 P α=1 for P . The independence test is converted to testing whether C P →A / ∈P B / ∈Q is a constant channel. Let ρ α = C P →A / ∈P B / ∈Q (ψ α ). We use SWAP tests to estimate the distances ρ 1 − ρ α 2 2 , and decide C P →A / ∈P B / ∈Q to be a constant channel if all distances are small.
According to Lemma 1, the SWAP test produces an estimate within error ε with probability no less than 1 − κ. Since Function 5 invokes the independence test O(n 2c+1 ) times, and each independence test invokes the SWAP test (Function 1) d 2 P = O(d 2c A ) times, there are O(n 2c+1 d 2c A ) number of SWAP tests. Therefore, the probability that all SWAP tests are within error ε is no less than 1 − Kn 2c+1 d 2c A κ for some constant K. Let κ = κ 0 /Kn 2c+1 d 2c A , and from now we assume all SWAP tests are within error ε, which has probability no less than 1 − κ 0 .
Let C (n) be the initial value of C in Function 5, and let C (k) , which has k input wires and k output wires, be the value of C in the (n − k + 1)-th recursive call of Function 5.
According to Function 1, proportional to a SIC-POVM for every subset P , we have λ −1 min = max P (d P + 1)/d P = O(1) [43]. The sample complexity is then and the computational complexity is O(T sample nc log d A ) from the realization of SWAP tests. Last, we show the bound of rank(C (k) ). Let A (k) and B (k) be the set of input and output wires C (k) , respectively. By the correctness Function 5, wires belonging to C (k) are before the wires not belonging to C (k) in the causal unravelling of C, and thus C ∈ Comb[(A (k) , B (k) ), (A (n) \A (k) , B (n) \B (k) )]. By Corollary 1, one has rank(C (k) ) ≤ rank(C) Bj ∈B (n) \B (k) d Bj / Ai∈A (n) \A (k) d Ai = rank(C (n) )d