Given black-box access to the input and output systems, we develop the first efficient quantum causal order discovery algorithm with polynomial query complexity with respect to the number of systems. We model the causal order with quantum combs, and our algorithm outputs the order of inputs and outputs that the given process is compatible with. Our algorithm searches for the last input and the last output in the causal order, removes them, and iteratively repeats the above procedure until we get the order of all inputs and outputs. Our method guarantees a polynomial running time for quantum combs with a low Kraus rank, namely processes with low noise and little information loss. For special cases where the causal order can be inferred from local observations, we also propose algorithms that have lower query complexity and only require local state preparation and local measurements. Our algorithms will provide efficient ways to detect and optimize available transmission paths in quantum communication networks, as well as methods to verify quantum circuits and to discover the latent structure of multipartite quantum systems.