Discretized quantum adiabatic process for free fermions and comparison with the imaginary-time evolution

Motivated by recent progress of quantum technologies, we study a discretized quantum adiabatic process for a one-dimensional free fermion system described by a variational wave function, i.e., a parametrized quantum circuit. The wave function is composed of $M$ layers of two elementary sets of time-evolution operators, each set being decomposed into commutable local operators. The evolution time of each time-evolution operator is treated as a variational parameter so as to minimize the expectation value of the energy. We show that the exact ground state is reached by applying the layers of time-evolution operators as many as a quarter of the system size. This is the minimum number $M_B$ of layers set by the limit of speed, i.e., the Lieb-Robinson bound, for propagating quantum entanglement via the local time-evolution operators. Quantities such as the energy $E$ and the entanglement entropy $S$ of the optimized variational wave function with $M<M_B$ are independent of the system size $L$ but fall into some universal functions of $M$. The development of the entanglement in these ansatz is further manifested in the progressive propagation of single-particle orbitals in the variational wave function. We also find that the optimized variational parameters show a systematic structure that provides the optimum scheduling function in the quantum adiabatic process. We also investigate the imaginary-time evolution of this variational wave function, where the causality relation is absent due to the non-unitarity of the imaginary-time evolution operators, thus the norm of the wave function being no longer conserved. We find that the convergence to the exact ground state is exponentially fast, despite that the system is at the critical point, suggesting that implementation of the non-unitary imaginary-time evolution in a quantum circuit is highly promising to further shallow the circuit depth.


I. INTRODUCTION
Currently realized and near-future expected quantum computing devices, called noisy intermediate-scale quantum (NISQ) devices [1], suffer various noise due to the poor gate fidelity and short coherent time so that the number of quantum gates as well as qubits available reliably in quantum devices is severely limited. Therefore, it is highly desirable to find quantum algorithms working efficiently on such a limited condition. One of the great challenges in quantum computing is to demonstrate quantum supremacy for practical calculations in quantum devices that can outperform the classical counterparts [1][2][3].
Quantum simulations of quantum many-body systems, such as the Hubbard model and quantum chemistry sys-tems, have been anticipated to be the most promising application for quantum computers [4]. One of the prominent algorithms specially in the NISQ era is the variational quantum eigensolver (VQE) [5,6], a quantumclassical hybrid algorithm, in which a variational wave function describing a quantum state is represented by a quantum circuit composed of parametrized quantum gates. In the VQE, the energy and often the derivatives of the energy with respect to the variational parameters are estimated on quantum computers and these quantities are used to optimize the variational wave function by minimizing the variational energy on classical computers.
In this regard, it has been pointed out that there occurs the barren plateau phenomena as potentially a serious issue in the VQE method [7]: if a random circuit is used, the derivative of a VQE wave function with respect to the parameters vanish as the number of qubits as well arXiv:2008.07168v1 [quant-ph] 17 Aug 2020 as quantum gates increases. Although this might not necessarily occur in general for any parametrized circuit, it is rather preferable to find a suitable circuit structure and reduce the number of parametrized gates necessary for representing a particular quantum state.
There have been several schemes proposed to improve and go beyond the plain VQE algorithm for quantum simulations of quantum many-body systems on the NISQ devices. One of the strategies is to systematically reduce the number of variational parameters in a parametrized quantum circuit with keeping the accuracy of the variational wave function. For example, the adaptive derivative-assembled pseudo-Trotter ansatz variational quantum eigensolver (ADAPT-VQE) [8] employs the unitary coupled cluster ansatz with generalized single and double excitations from a single Slater determinant reference, but the parametrized quantum gates are additively selected, one at each iteration, in the collection of one-and two-body operator pool by searching the appropriate gate that gives the largest gradient, and hence the optimization is performed only for the selectively accumulated quantum gates. In a symmetryadapted VQE scheme, the symmetry of Hamiltonian is imposed in the variational wave function to reduced the number of parametrized gates at the expense of introducing a non-unitary projection operator that is treated partly as postprocessing on classical computers [9]. A non-orthogonal VQE scheme is a multireference version of the VQE algorithm where a generalized eigenvalue problem in a subspace spanned by a collection of the parametrized variational wave functions is solved, in addition to optimizing the variational parameters [10]. A similar idea of expanding a quantum subspace within the VQE scheme is also proposed in Refs. [11][12][13][14].
Another strategy is to employ the imaginary-time evolution that is non-unitary. Recently, Motta et al. have proposed a quantum imaginary-time evolution (QITE) algorithm [15], in which a local infinitesimally small imaginary-time evolution operator, including the normalization factor of the imaginary-time evolved quantum state, is mapped to a non-local unitary real-time evolution operator by solving a linear system of equations on classical computers to properly parametrize the non-local unitary operator that approximately reproduces the local non-unitary imaginary-time evolution operator [15][16][17][18]. Note that the parameters in a non-local unitary operator here are determined by solving a linear system of equations, not by optimizing a cost function as in the VQE scheme. It is also interesting to note that in Ref. [15] they have used the QITE algorithm to expand a quantum subspace by constructing non-orthogonal Krylovsubspace basis states and proposed a quantum version of a Lanczos-like algorithm. In this regard, an imaginarytime evolution is not necessarily required to generate a Krylov subspace, as demonstrated in Refs. [19,20] by using a real-time evolution. Very recently, a quantum version of the power method is proposed to generate a Krylov subspace [21].
Considering application for the NISQ devices, it is crucial to find a way, based on some guiding principle, of designing a quantum circuit ansatz, ideally with linear depth or less, that can efficiently represent a quantum state of interest. It is easily shown mathematically that the imaginary-time evolution can yield the exact solution of a ground state in the limit of long-time evolution, provided that the imaginary-time evolution is treated exactly. The Lanczos method also guarantees to converge to a ground state with a desired accuracy as the dimension of a Krylov subspace is increased. Although these methods are well established classically, their quantum versions are still under development, as described above. There have been various circuit ansatze proposed in the VQE scheme such as a unitary coupled cluster ansatz for mostly quantum chemistry application [5,6,22,23] and a hardware efficient ansatz [24]. These ansatze can represent any quantum state in principle by increasing the number of gates and thus the circuit depth [25], and have been implemented in the NISQ devices with success particularly for small molecules [24,[26][27][28].
Here, in this paper, we shall focus on a circuit ansatz realized by discretizing a quantum adiabatic process from an initial product state to a final state corresponding to a ground state of a Hamiltonian to be solved [29][30][31][32]. This is inspired by the quantum approximate optimization algorithm (QAOA) for combinatorial optimization problems that are represented as an Ising model [33]. In this paper, this circuit ansatz is called a discretized quantum adiabatic process (DQAP) ansatz. An advantage of a DQAP ansatz is the fact that a circuit constructed by the DQAP can yield the exact ground state without any parametrization in the continuous circuit limit because of the quantum adiabatic theorem [34][35][36][37]. However, the convergence of a DQAP ansatz with the finite number of discretization steps is unknown in general and this is the main issue addressed in this paper.
We thereby study a DQAP ansatz for free fermions on a one-dimensional lattice at half filling. This is an ideal system to analyze a DQAP ansatz because a quantum state evolved by the DQAP with an initial state described by a single Slater determinant state can still be described by a single Slater determinant state, and therefore one can keep track of each occupied single-particle orbital in the Slater determinant state during the DQAP. In the DQAP ansatz considered here, we first prepare as the initial state a product state of local bonding states formed on neighboring sites, which can be described by a single Slater determinant state. We then let the state evolve forward via the DQAP by repeatedly applying layers of two elementary sets of local time-evolution operators (see Fig. 1), where the evolution time in each layer of the time-evolution operators is treated as a variational parameter so as to minimize the variational energy. We examine how the state described by the DQAP ansatz is evolved with increasing the number M of layers by monitoring the variational energy, single-particle orbitals in the Slater determinant state, the entanglement entropy, and the mutual information.
We find that the exact ground state is attained by applying the layers of time-evolution operators as many as a quarter of the system size, which is the minimum number M B of layers necessary to entangle the entire system by the local time-evolution operators, corresponding to the Lieb-Robinson bound for the propagation of quantum information [38]. In contrast, the DQAP ansatz with the number M of layers less than M B , thus not describing the exact ground state, represents another series of quantum states in that physical quantities such as the energy and the entanglement entropy evaluated for these states are independent of the system size but scale with M , which indicates that the entanglement carried by the DQAP ansatz with a finite number of layers is bounded, as in the case of the matrix product states with a finite bond dimension [39,40]. We also find that the optimized variational parameters in the DQAP ansatz converge systematically to a smooth function of the discretized time, which thus provides the optimized scheduling function for the quantum adiabatic process. Furthermore, we show that the effective total evolution time of the optimized DQAP ansatz with M B layers of the local time-evolution operators is proportional to the system size. For comparison, we also investigate the imaginary-time evolution of the DQAP ansatz, which can still be described by a single Slater determinant state. We find that the convergence to the ground state is exponentially fast with respect to the number of layers of the local imaginary-time evolution operators, despite that the system is at the critical point where the one-particle density matrix decays algebraically with distance.
The rest of this paper is organized as follows. We first describe the free fermion model and establish the notation used throughout this paper in Sec. II A, and introduce the DQAP ansatz in Sec. II B. We then provide the analytical formulae for various quantities such as the one-particle density matrix, entanglement entropy, and mutual information, and also explain the optimization method in Secs. II C-II E. The numerical results for the DQAP ansatz are given in Sec. III and these results are compared with those for the imaginary-time counterpart in Sec. IV. We then conclude the paper with brief discussion in Sec. V. The numerical details of the optimized parameters are discussed in Appendix A and details of the entanglement entropy with respect to the one-particle density matrix are explained in Appendix B. Throughout the paper, we set the reduced Planck constant = 1.

II. MODEL AND METHOD
In this section, we first define the free fermion model with matrix notation in Sec. II A, and introduce a discretized quantum adiabatic process to construct a variational ansatz in Sec. II B. We then summarize the analytical formulae for the variational ansatz in the free fermion case in Sec. II C. The optimization method to optimize the variational parameters is described in Sec. II D. To discuss the entanglement property, we also derive the analytical formulae of the reduced density matrix, entanglement entropy, and mutual information for a free fermion wave function in Sec. II E.

A. Model
The free fermion system considered in this paper is described by the following Hamiltonian: whereĉ x (ĉ † x ) denotes the annihilation (creation) operator of a fermion at site x ∈ {1, 2, · · · , L}. For convenience, we represent the Hamiltonian in Eq. (1) aŝ whereĉ † (ĉ) is a L dimensional row (column) vector of the fermion operators given bŷ and T is an L × L matrix whose elements are given by Let U be a unitary operator which diagonalizes the matrix T as where E is the L × L diagonal matrix whose diagonal elements are the eigenvalues of T : E = diag(E 1 , E 2 , · · · , E L ).
Here, we assume E n (n = 1, 2, · · · , L) in ascending order, i.e., E 1 ≤ E 2 ≤ · · · ≤ E L . The ground state of the Hamiltonian in Eq. (1) for N fermions is represented as where |0 is the vacuum of fermions and Ψ is an L × N matrix obtained by extracting the first N columns from U . [ĉ † Ψ] n indicates the nth element of the N dimensional row vector c † Ψ.
It is important to note that the row index x of [Ψ] xn indicates the site index while the column index n is the index labeling the single-particle state with the singleparticle energy E n obtained by diagonalizing the matrix T . [ĉ † Ψ] n in Eq. (5) thus corresponds to the nth single-particle orbital that composes the Slater determinant state of the ground state |ψ . Hereafter, we simply refer to the column vectors of Ψ as single-particle orbitals.
B. Variational ansatz based on a discretized quantum adiabatic process The quantum adiabatic process is a quantum process following the quantum adiabatic theorem [34][35][36][37], in which a slowly driving system in time from an initial eigenstate, the ground state of HamiltonianĤ i , stays in the instantaneous eigenstate of the time evolving Hamil-tonianĤ(τ ) at time τ and finally reaches to the ground state of HamiltonianĤ f . Here the time evolving Hamil-tonianĤ(τ ) is expressed aŝ with s i (τ ) and s f (τ ) being the scheduling functions that are smooth and satisfy the conditions: where τ i and τ f denote the initial and final times of the process, respectively. The final state |ψ(τ f ) at τ = τ f after the time evolution is thus given as whereÛ(τ, τ i ) is the time-evolution operator obtained by solving the Schrödinger's equation withÛ(τ i , τ i ) = 1 and |ψ i is the ground state ofĤ i . It is well known that for a sufficiently long time τ f − τ i , implying a slow driving dynamics, the initial state |ψ(τ i ) = |ψ i is adiabatically transformed into the ground state of the final HamiltonianĤ f through this adiabatic process if there is a finite energy gap between the ground state and the excited states ofĤ(τ ) for all τ [41,42]. A quantum adiabatic process is a real-time dynamics governed by the unitary time-evolution operator in Eq. (8). This should be contrasted with the case of the imaginary-time evolution where the imaginary-time evolution operator is no longer unitary. Since all operations in quantum computers are composed of unitary gate operations, a quantum adiabatic process would be a natural principle to follow in constructing a circuit ansatz for obtaining a ground state of a Hamiltonian in a quantum circuit.
We shall now consider the case where the final Hamil-tonianĤ f in Eq. (6) is composed of a set of termsV p (p = 1, 2, · · · , P ):Ĥ f = P p=1V p (10) such that in general [V p ,V p ] = 0 when p = p . Here, we assume that eachV p consists of a set of operatorŝ where all operatorsÔ (p) q commute with each other for given p: In addition, we assume that the initial HamiltonianĤ i in Eq. (6) is given by one ofV p 's inĤ f and here we consider Then, the time-evolution operator is written in the following form: where in order to reproduceÛ(τ f , τ i ) in the limit of M → ∞. This is the most naive discretization procedure of time in the time-evolution operatorÛ(τ f , τ i ) and the ground state ofĤ f is obtained by operatingÛ(τ f , τ i ) to the initial state |ψ i as in Eq. (8). The simplest scheduling functions s i (τ ) and s f (τ ) that satisfy the conditions given in Eq. (7) is a linear scheduling, i.e., where T = τ f − τ i [42,43]. In this case, Inspired by the quantum adiabatic process described above, here we instead consider, as a variational ansatz for the ground state ofĤ f , the following state with a finite value of M : where θ = {θ m } M m=1 are assumed to be variational parameters determined by minimizing the variational energy. This variational state |ψ M (θ) in Eq. (21) is referred to as a discretized quantum adiabatic process (DQAP) ansatz.
Let us now illustrate the DQAP ansatz for the free fermion system given in Eq. (1). For simplicity, we assume that the system is one dimensional and the final HamiltonianĤ f is given bŷ where γ sets the boundary conditions: γ = 1 for the periodic boundary conditions (PBC) and γ = −1 for the anti-periodic boundary conditions (APBC). We also assume that the number L of sites is even and the number N of fermions is at half filling, i.e, N = L/2. In what follows, we set t = 1 as a unit of the energy.
FIG. 1. Schematic representation of the discretized quantum adiabatic process (DQAP) ansatz |ψM (θ) , defined in Eq. (27), for the one-dimensional free fermion system. Horizontal black solid lines indicate the local quantum states at sites x (i.e., qubits). Green blocks denote the local timeevolution operators of the form exp[iθt(ĉ † xĉx +ĉ † x ĉx)]. The boundary terms (blocks connecting the top and bottom black solid lines) involve generally non-local operations in the qubit representation, but they can be eliminated for the special cases discussed in the text. The initial state |ψi is given in Eq. (25). Light yellow squares indicate the local bonding The schematic representation of this state is shown in Fig. 1.
We shall now discuss how this DQAP ansatz can be described in the qubit representation for quantum computing. First of all, a fermion system can always be mapped in the qubit representation through, e.g, the Jordan-Wigner transformation [44]: σ ± x = (X x ±iŶ x )/2, and {X x ,Ŷ x ,Ẑ x } are the Pauli operators (i.e., gates) at qubit x. Notice that [σ ± x ,K ( †) (x)] = 0 and from Eq. (28)ĉ † xĉx = 1 2 (Ẑ x + 1). With this transformation, any local fermion operator acting up to nearestneighbor sites in a one-dimensional system can be represented by Pauli operators without introducing the sign factors due to the Jordan-Wigner stringK(x), suggesting that the fermion representation is trivially equivalent to the qubit representation, except for the boundary terms. Indeed, the sign factors at the boundary are also cancelled if APBC (PBC) is imposed when N is even (odd).
In this case, the one-dimensional free fermion system in Eq. (22) under both PBC and APBC is mapped onto the spin-1/2 XY model withσ ± L+1 =σ ± 1 (i.e, PBC). Here,σ + x (σ − x ) represents the local operator to flip the qubit state from |1 x (|0 x ) to |0 x (|1 x ), but not the other way around, where |σ x (σ = 0, 1) denotes the local state at qubit x in the Pauli z basis. As shown in Fig. 2(a), the local time-evolution op- can be implemented in a quantum circuit [45][46][47]. It should also be noted that the condition of N being even (odd) for APBC (PBC) corresponds to the closed shell condition in the free fermionic system.
With the Jordan-Wigner transformation, the fermion vacuum state |0 is mapped to L x=1 |1 x , and therefore the initial state |ψ i in Eq. (25) can be mapped in the qubit representation onto i.e., a product state of spin-triplet states. As shown in Fig. 2 Here,Ĉ x (X x ) denotes the control-NOT gate acting on qubit x with the control qubit at qubit x, andĤ x indicates the Hadamard gate acting on qubit x. Finally, we briefly note that for a more general fermion system in higher spatial dimensions with a long-range hopping, the phase factors due to the Jordan-Wigner strings cannot be canceled and yield many-body interactions in the qubit representation. In principle, these many-body interactions can be treated as two-qubit operations by using, for example, the perturbative gadgets [48]. However, these techniques introduce additional sources of errors. Therefore, we leave the general cases for a future study and focus here on the one-dimensional system.

C. Useful properties of the DQAP ansatz for free fermions
The DQAP ansatz for the free fermion system introduced in the previous section has generally the following form: whereŴ k (k = 1, 2, · · · , K) is a Hermitian single-particle operator given byŴ and |ψ 0 is a ground state of an N -fermion system defined by the following single-particle Hamiltonian: thus representing a Slater determinant state of N fermions. We can now easily show that |ψ(θ) in Eq. (33) is more compactly written as where and Ψ 0 is an L × N matrix such that the nth column of Ψ 0 is the eigenstate of W 0 with the nth lowest eigenvalue. The derivation of Eq. (36) can be found in Ref. [49]. Equations (36) and (37) imply that a state initially prepared as a single Slater determinant state evolves in time, realized by repeatedly applying the unitary timeevolution operators, into a state that can still be represented as a single Slater determinant state. Therefore, we can even discuss the time evolution of each constituent single-particle orbital in the Slater determinant state. It is also ready shown that the overlap between two N -fermion states |ψ and |φ is calculated as where |ψ is an N -fermion state given in Eq. (5) but for any Ψ and with Φ being an L × N matrix. We can also show the following useful formula: where tr[A] indicates the trace of a matrix A and δ xx is an L×L matrix whose elements are given by [δ xx ] x1x2 = δ xx1 δ x x2 . We can furthermore derive that, for example, which is simply the Wick's theorem.

D. Optimization method
In this paper, we employ the natural gradient method to optimize the variational parameters in the variational wave function. Here we briefly summarize this optimization method for the DQAP ansatz.
The natural gradient method was originally introduced in the context of machine learning [50,51]. However, essentially the same method has also been independently proposed to optimize a many-body variational wave function [52] and it has been successfully applied to various systems in quantum chemistry and condensed matter physics [53,54]. This method has also been proposed recently in the context of quantum computing as a way to optimize a parametrized quantum circuit [55] and is nicely summarized in Ref. [56].
It is now well known that there are several ways to derive this optimization method [57]. A simple way is the formulation based on an infinitesimal imaginary-time evolution in the variational parameter space. In this case, we determine the new variational parameters θ new = θ + δθ so as to satisfy where |ψ(θ) is given in Eq. (33) with K variational parameters θ = {θ 1 , θ 2 , · · · , θ K },Ĥ is the Hamiltonian to be solved, in our case, given in Eq. (1), and δβ is an infinitesimal real number. Assuming that the variational parameters are all real, δθ is then determined as where d(|ψ , |φ ) is a distance between two quantum states |ψ and |φ and is given by i.e., essentially the same as the fidelity, assuming that the two states |ψ and |φ are not generally normalized. Expanding d 2 (|ψ(θ + δθ) , (1 − δβĤ)|ψ(θ) ) up to the second order of δθ and δβ, we obtain the following quadratic form: where δθ in the right hand side is a K-dimensional column vector with the kth element being δθ k , S is a K ×K matrix given by with f is a K-dimensional column vector given by andĒ 2 is the variance of the HamiltonianĤ given bȳ Here, we assume that |ψ(θ) is not normalized and thus these formulae can be used in general cases. Note also that S is Hermitian, i.e., S † = S. The stationary point of the quadratic equation given in Eq. (45) is now easily obtained by solving the following linear equation: Notice that since S + S * and f + f * are both real, the solution δθ t = (δθ 1 , δθ 2 , · · · , δθ K ) is guaranteed to be also real. We can thereby obtain the new variational parameters θ new = θ + δθ by solving the above linear equation, in which δβ is learning rate and can be chosen properly.
We can now easily show that Since S + S * = 2Re[S] is a positive semidefinite matrix [9], δE ≤ 0 as long as δβ > 0. We should also note that if we expand the following quantity: the matrix S defined in Eq. (46) naturally appears. Therefore, S can be regarded as a metric tensor for the distance d(|ψ , |φ ) in the parameter space θ.
Using Eqs. (38) and (40), we can derive explicitly the forms of S and f , respectively, for the variational state given in Eq. (33) as and where we have used that Ψ † K Ψ K = I N and I N is the N -dimensional unit matrix. This condition corresponds to the fact that the single-particle orbitals in Ψ K are orthonormalized, and the generalization to the case where they are not orthonormalized is described in Sec. IV. ∂ k Ψ K is an L×N matrix defined as the first derivative of Ψ K given in Eq. (37) with respect to the kth variational parameter θ k , i.e., Finally, notice that the update formula given in Eq. (50) can be regarded as an extension of the steepest descent algorithm that corresponds to the case when the metric tensor S is the unit matrix, indicating that the optimization method described here cannot exceed the limitation of the locality of the search space in general. However, we find that this is not a problem in our case since we can easily obtain the optimal results without any difficulty. The detail of this point is found in Appendix A.

E. Entanglement entropy for free fermions
The entanglement von Neumann entropy is a measure to quantify the quantum entanglement between a subspace and its complement of a quantum state, and has been used to characterize various quantum states. The formula is quite simplified for the free fermion systems and here we briefly outline how the entanglement entropy is calculated.
Let A = {x 1 , x 2 , · · · , x L A } be a subset of sites (the number of sites in A being L A ) that are picked out of the all sites U = {1, 2, · · · , x, · · · , L}. Let B be the complementary subspace of A: B = A. We also assume that |ψ is a normalized quantum state and can be represented by using the basis on U. The reduced density matrixρ A of the subspace A is given bŷ where Tr B indicates the trace over all basis defined on the subspace B. The entanglement entropy S A of the subspace A is defined by using this reduced density matrixρ A as where Tr A is the trace over all basis defined on the subspace A.
Notice first that the expectation value of any physical quantityÔ A defined on the subspace A can be obtained by using the reduced density matrixρ A as For the fermion system,Ô A is generally composed of a product ofĉ x andĉ † x with x ∈ A. Therefore, when |ψ is a single-particle state, we can use the Wick's theorem [see, for example, Eq. (41)]. This implies thatρ A can be written asρ whereĉ † A andĉ A are similar to those in Eq. (3) but the elements here are fermion operatorsĉ † x andĉ x with x ∈ A [58]. Indeed, one can derive the matrix Γ directly from a given single Slater determinant state |ψ [58,59]. Here, we shall follow a different route [60].
Since the Wick's theorem can decompose the expectation value of any operator into a product of one-particle density matrices, we can determine Γ by equating the expectation values of the single-particle operatorĉ † xĉx , i.e., To this end, let us introduce the following L A × L A one-particle density matrix: whereĉ * A (ĉ t A ) is the matrix transpose ofĉ † A (ĉ A ) and |ψ is a single Slater determinant state given by the form of Eq. (5). Using Eq. (40), each element of D A can be obtained as where V denotes the unitary matrix composed of the eigenvectors of matrix D A and ∆ is the diagonal matrix whose diagonal elements correspond to the eigenvalues δ l (l = 1, 2, · · · , L A ) of matrix D A : ∆ = diag(δ 1 , δ 2 , · · · , δ L A ). Let us also define the following L A × L A matrix: assuming thatρ A is given in Eq. (59). We then obtain that where I L A is the L A -dimensional unit matrix, U A is the unitary matrix composed of the eigenvectors of matrix Γ, and Λ is the diagonal matrix whose diagonal elements correspond to the eigenvalues λ l (l = 1, 2, · · · , L A ) of matrix Γ: Λ = diag(λ 1 , λ 2 , · · · , λ L A ). Because of Eq. (58), we now impose that D A = D A . Comparing Eqs. (62) and (64), we obtain that and Therefore, we finally find that Giving the form of the reduced density matrixρ A in Eq. (59), the entanglement entropy S A defined in Eq. (57) can now be written as Using Eq. (66), we can show that the entanglement entropy S A is simply given as Therefore, the entanglement entropy S A is determined solely from the eigenvalues of the one-particle density matrix D A . Note that the eigenvalues of D A are bounded as 0 ≤ δ l ≤ 1. In the context of quantum chemistry, the eigenvectors of D A are called the natural orbitals and the eigenvalue δ l corresponds to the density of each natural orbital. Since the L A × L A matrices inside the trace in Eq. (69) are diagonal, we can discuss separately the individual contribution of the natural orbitals to the entanglement entropy S A . For example, the contribution to S A is maximum when δ l = 0.5, while it is minimum when δ l = 0 or δ l = 1. This implies that when δ l = 0.5, the corresponding natural orbital in the subsspace A are highly hybridized with orbitals in the subspace B, giving an intuition of the quantum entanglement in the free fermion system. As described above, the entanglement entropy S A is a measure to quantify the quantum entanglement between subspaces A and B = A. Instead, it is often required to discuss how the quantum state is entangled between a subspace A ⊂ U and another subspace B ⊂ U with A ∩ B = ∅ and A ∪ B = U. One of the quantities for this purpose is the mutual information I A,B defined by We consider a special case when A = {x} and B = {x }, and the mutual information I x,x for this special case is There are several remarks on I x,x . First of all, D {x} = D x = ψ|ĉ † xĉx |ψ , which is the density of fermions at site x. Therefore, if the system is homogenous, D x is independent of x and D x = N/L. When the system is at half filling, D x = 0.5 and thus S {x} = ln 2, which is the maximum value of the entanglement entropy for a single site. Second, I x,x is determined by S {x,x } , which can be calculated form the eigenvalues of the one-particle density matrix Since the diagonal term is 0.5 when the system is homogenous at half filling, the off-diagonal elements determine the value of I x,x . For example, if | ψ|ĉ † xĉx |ψ | = 0.5, we find that S {x,x } = 0 and thus I x,x = 2 ln 2. In contrast, if ψ|ĉ † xĉx |ψ = 0, we find that S {x,x } = 2 ln 2 and thus I x,x = 0.

III. NUMERICAL RESULTS
Here, we show the results of numerical simulations for the one-dimensional free fermion system described in Eq. (22) and examine how the DAQP ansatz |ψ M (θ) given in Eq. (27) approaches to the exact ground state with increasing the number M of layers of the local timeevolution operators (see Fig. 1). We focus on the fermion density at half filling, i.e., N = L/2 and use the natural gradient method described in Sec. II D to optimize the variational parameters θ in the DAQP ansatz |ψ M (θ) .

A. Convergence of ground state energy
We optimize the variational parameters θ = {θ } in the DAQP ansatz |ψ M (θ) given in Eq. (27) so as to minimize the variational energy for a given system size L. In order to check the convergence of the variational energy, Fig. 3(a) shows the energy difference ∆E = E M (L) − E exact (L) from the exact energy E exact (L) for various system sizes L as a function of M . Here, we use L = 4n L (n L : integer) with APBC and thus the closed shell condition is satisfied for the ground state (see Sec. II B). As shown in Fig. 3(a), ∆E monotonically decreases with increasing M and we obtain that ∆E = 0 within the machine precision exactly at M = L/4 for all values of L studied [61].
To better understand this observation, let us examine closely how the expectation value of the energy for the DQAP ansatz |ψ M (θ) in Eq. (73) is evaluated. For this purpose, we should notice that the energy expectation value is essentially given simply by the sum of terms ψ M (θ)|ĉ † xĉx+1 |ψ M (θ) (and also ψ M (θ)|ĉ † x+1ĉ x |ψ M (θ) but it is basically the same as ψ M (θ)|ĉ † xĉx+1 |ψ M (θ) for the purpose of the discussion here) over all x's. Therefore, it is adequate to consider each term separately. Because of the form of construction for the DQAP ansatz |ψ M (θ) , there are two different cases of ψ M (θ)|ĉ † xĉx+1 |ψ M (θ) : the operatorĉ † xĉx+1 acts (i) over two neighboring local time-evolution opera- , as schematically shown in Fig. 4(a), and acts (ii) only on a single local time-evolution operator e itθ (M ) , as shown in Fig. 4(b).
Let us first consider case (i). In this case, wherê assuming that L ≥ 4M + 2. Namely, one can eliminate many of the local unitary time-evolution operators in the expectation value due to the cancellations of the left and right sides of the product, as illustrated in Fig. 4(a). The number i of sites (i.e., qubits) that contribute to the local expectation value ψ M (θ)|ĉ † xĉx+1 |ψ M (θ) is linearly dependent on M : i = 4M + 4 [62]. This implies that the propagation of quantum entanglement via the local timeevolution operators is bounded in space and this boundary forms a causality-cone like structure shown schematically in Fig. 4(a). This upper limit on the propagation speed is known as the Lieb-Robinson bound [38].
In the case of case (ii), we can also evaluate the local expectation value ψ M (θ)|ĉ † xĉx+1 |ψ M (θ) in the same manner as in Eq. (74) except that noŵ andV assuming that L ≥ 4M . Therefore, in this case, the number ii of sites that contribute to the local expectation value ψ M (θ)|ĉ † xĉx+1 |ψ M (θ) is also linearly dependent on M : ii = 4M + 2 [62]. This sets the boundaries of a causality-cone like structure in Fig. 4(b), within which the quantum entanglement is developed.
In order to reach the exact ground state energy for a given system size L, i and ii have to be equal to or exceed the system size L, which corresponds to with z being the smallest integer greater than or equal to z. This condition is independent of the boundary conditions because Eqs. The same conclusion is also reached in the case when the PBC is employed [61]. We should also note that as indicated in Fig. 4, the causality-cone like structure of the local unitary timeevolution operators contributing to the local expectation value does not depend on the system size L. As a consequence, we expect that the optimized variational energy per site would not depend on L as along as L > 4M +2 (= ii ). Indeed, as shown in Fig. 3(b), the optimized variational energies per site with a given value of M are exactly the same for different values of L until M reaches to the boundary at M = L/4 for APBC [63], where the variational energy abruptly changes to the exact value for the system size L. Moreover, we find that the optimized variational energy per site for M < (L − 2)/4 under PBC is identical to that for the same M (but M < L/4) under APBC. We should note that a similar analysis for the transverse-field Ising model has also been reported in Ref. [30].

B. Time-evolution of single-particle orbitals
We now explore how the DQAP ansatz |ψ M (θ) is evolved by applying the local time-evolution operators. Following the argument in Sec. II C, the DQAP ansatz |ψ M (θ) given in Eq. (27) can be written as where Here, V 1 and V 2 are L × L matrices representingV 1 and V 2 given in Eqs. (23) and (24), respectively, i.e., with and and Ψ i is the L × N matrix representing the initial state |ψ i in Eq. (25), i.e., with It should be noted that since V 1 and V 2 are both block diagonal matrices with each block being a 2 × 2 matrix, these can easily be exponentiated as and It is also apparent from Eq. (86) that each column vector in Ψ i corresponds to a single-particle orbital, representing the local bonding state 1 √ 2 (ĉ † 2x−1 +ĉ † 2x )|0 in this case given in Eq. (26), which constitutes the Slater determinant state |ψ i for N fermions. Therefore, we can now clearly understand that the time-evolved state |ψ M (θ) from a state initially prepared as a single Slater determinant state |ψ i can still be represented as a single Slater determinant state, in which each single-particle orbital is given by each column vector of Ψ M . We can thus examine the time evolution of each single-particle orbital in the Slater determinant state, which is described by Eqs. (81), (87), and (88).
For this purpose, we first introduce the following L×N matrix: for m = 0, 1, 2, · · · , M , where the variational parameters } are determined so as to minimize the variational energy for |ψ M (θ) , and thus Φ M = Ψ M . We also set that Φ 0 = Ψ i . There are two elemental properties of Φ m . First, the single-particle orbitals in Φ m are mutually orthonormalized. This is simply because of the consequence of the unitary evolution: Second, it is apparent by construction in Eq. (89) that, apart from the phase factor due to the boundary conditions (i.e., in the case of APBC), a single-particle orbital in Φ m is transformed to other single-particle orbitals by the translation of two lattice space, i.e., [Φ m ] x,n = [Φ m ] x±2,n±1 , where n = 0 and N + 1 correspond to N and 1, respectively. Therefore, the single-particle orbitals in Φ m are associated with the Wannier orbitals with a unit cell of two lattice space.
Let us now introduce the spatial extent d m (in unit of lattice constant) of a single-particle orbital in Φ m , i.e., d m being the number of consecutive nonzero elements in each column of Φ m . It is obvious form Eq. (86) that d 0 = 2 for Φ 0 = Ψ i . Without knowing the explicit values of the variational parameters θ, we can readily show that the spatial extent of a single-particle orbital increases by four each time applying matrices e −iθ (m) 1 V1 and e −iθ (m) 2 V2 given in Eqs. (87) and (88), respectively, i.e., d m = d m−1 + 4. Therefore, the spatial extent of a singleparticle orbital in Φ m is generally given as d m = 4m + 2 for our initial matrix Φ 0 = Ψ i . Consequently, the spatial extent d m of a single-particle orbital in Φ m exceeds (reaches) the system size L at m = L/4 [m = (L − 2)/4] for APBC (PBC) where we choose L = 4n L (L = 4n L +2) with n L integer. In other words, in order for the singleparticle orbitals in Ψ M = Φ M to extend over the entire system, the smallest number M of layers in Ψ M is L/4 [(L − 2)/4] for APBC (PBC), which is in good accordance with the results in Fig. 3 and the discussion in Sec. III A. This is understood because the spatial extent d m of the single-particle orbitals essentially sets the limit of the propagation of quantum entanglement in the DQAP state. Figure 5 shows the numerical results of the time evolution of a single-particle orbital in Ψ M of the DQAP ansatz |ψ M (θ) , for which the variational parameters θ are optimized for L = 90 with M = (L−2)/4 under PBC, thus representing the exact ground state. Initially, the single-particle orbital is spatially localized at sites x = 1 and 2, and propagates gradually in time (i.e., increasing m) by splitting a wave into the two opposite directions, finally reaching each other at m = M when the spatial   extent d m of the single-particle orbital becomes as large as the system size L.
We should note here that the single-particle orbitals in Ψ M are not uniquely determined. Instead, an L × L matrix Ξ M given by is invariant for all sets of single-particle orbitals which represent the exact ground state. The expectation value of any physical operator evaluated for |ψ M (θ) is the same, despite that Ψ M is not uniquely determined, as long as Ξ M is the same for different Ψ M . This can be easily proved from Eq. (40) because the single-particle orbitals are orthonormalized, i.e., Ψ † M Ψ M = I N . It is also apparent that Ξ M is invariant under the transformation where Q is an N × N unitary matrix. Starting with dif-ferent initial variational parameters, the numerical optimization of the variational parameters in |ψ M (θ) might provide different sets of optimized variational parameters and thus different Ψ M 's. We indeed find several sets of single-particle orbitals with different single-particle orbital shapes, which nonetheless constitute the exact ground state, and all of them give the same value of Ξ M . However, we note that all of these sets of single-particle orbitals are time evolved as those shown in Fig. 5, and they extend over the entire system at m = M = (L−2)/4 for PBC.
We shall now consider the number of independent matrix elements in an L × N complex matrix Ψ when the exact ground state is constructed in the form |ψ = N n=1 [ĉ † Ψ] n |ψ i where |ψ i is given in Eq. (25). To be specific, we assume that L = 4n L + 2 (n L : integer) with PBC at half filling, i.e., N = L/2. We should first recall that [Ψ] x,n represents the nth single-particle orbital at site x. Because a single-particle orbital can be mapped to other single-particle orbital by the translation of two lattice space, i.e., [Ψ] x,n = [Ψ] x±2,n±1 , there are L independent complex elements in Ψ. In addition, there exists the reflection symmetry at the center of bond, i.e., [Ψ] x,n = [Ψ] −x+4n−1,n , which reduces the number of independent complex elements in Ψ down to L/2. Furthermore, the orthonormality of the single-particle orbitals, i.e., Ψ † Ψ = I N , yields n L + 1 independent equations and thus there are 3n L + 1 independent real elements in Ψ.
Next, we shall consider the transformation of Ψ by Q, i.e., Ψ → Ψ = ΨQ, as dicssused above in Eq. (92). Assuming that Ψ has the same translational and reflection symmetries as in Ψ, we can show that the matrix elements of Q are also related to each other, similarly to the matrix elements of Ψ. Thus, the independent complex matrix elements in Q is n L + 1. In addition, the unitarity of Q yields n L + 1 independent equations and therefore there are n L + 1 independent real elements in Q. This suggests that, among 3n L + 1 independent real elements in Ψ, n L + 1 elements are redundant. Therefore, there are 2n L = (L − 2)/2 independent real elements that represent Ψ. It is interesting to notice that this number coincides with the number of the variational parameters θ = {θ In the case of APBC with L = 4n L , exactly the same argument follows except that now there are n L independent equations generated due to the orthonormality of the single-particle orbitals in Ψ and the unitary matrix Q contains n L independent real elements. Therefore, there are 2n L = L/2 independent real elements in Ψ. This number also coincides with the minimum number of variational parameters θ in the DQAP ansatz |ψ M (θ) with M = L/4 that can represent the exact ground state.

C. Entanglement entropy
Next, we shall examine the entanglement property of the DQAP ansatz |ψ M (θ) . Figure 6 shows the entanglement entropy S A of the optimized DQAP ansatz |ψ M (θ) as a function of M for several different system sizes L under APBC. Here, the variational parameters θ are optimized for each M , and the bipartition is assumed to be the half of the system, namely, with the size L A of subsystem A being L/2. For the bipartitioning, we consider only the case where the system is divided into the two subsystems by not breaking any local bonding state in the initial state |ψ i , as shown schematically in Fig. 7. As shown in Fig. 6, we find that the entanglement entropy S A for 4M ≤ L A is independent of the system size L and falls into a smooth "universal" curve of M . On the other hand, the entanglement entropy S A starts to deviate from this universal curve for 4M > L A and approaches to the exact value at M = L/4 for APBC [M = (L − 2)/4 for PBC]. These features are captured schematically in Fig. 7. The partitioning effect can propagate via the local unitary time-evolution operators into the inside of subsystem A up to 2M lattice space (taking also into account the entanglement of a local bonding state in |ψ i ) from each boundary of the partitioning, and thus this maximum propagation limit forms a causality-cone line structure centered at each partitioning boundary (see Fig. 7). The two causality cones cross each other when 4M > L A , and this is when the entanglement entropy S A deviates from the universal curve found in Fig. 6.
Let us discuss the results for 4M ≤ L A . In this case, we find that the entanglement entropy S A is exactly the same for different system sizes L and thus different sizes L A of subsystem A. This implies that the entanglement entropy S A is independent of the size L A of subsystem A, as long as the partitioning boundaries are separated long enough. In other words, the entanglement carried by the DQAP ansatz |ψ M (θ) with a finite M is bounded, as in the matrix product states with a finite bond dimension [64,65].
So far, we have assumed that the size LĀ of the com-plementĀ of subsystem A is the same as the size L A of subsystem A. However, we should note that the results of the entanglement entropy S A for 4M ≤ L A shown in Fig. 6 remain exactly the same even when we enlarge the size ofĀ, provided that L A ≤ LĀ. Thus, the entanglement entropy S A for 4M ≤ L A and L A ≤ LĀ is determined solely by the number M of layers in the DQAP ansatz |ψ M (θ) . Note also that, considering S A = SĀ, the smaller subsystem determines the value of M until which the entanglement entropy follows the universal curve.
We have performed similar calculations for the systems under PBC and found that, independently of the system size L, the entanglement entropy S A of the optimized DQAP ansatz |ψ M (θ) for the system under PBC is exactly the same as that of the optimized DQAP ansatz |ψ M (θ) for the system under APBC, provided that i) 4M ≤ L A , ii) L A ≤ LĀ, and iii) the partitioning boundaries do not break any local bonding state in the initial state |ψ i . Namely, the entanglement entropy S A falls into the universal curve of M shown in Fig. 6, independently of the system size and the boundary conditions, as long as the three conditions above are satisfied. This is similar to the observation of the optimized variational energy E M (L), where ∆ε = E M (L)/L − ε ∞ falls into the universal curve of M , as shown in Fig. 3(b), independently of the system size and the boundary conditions, as long as M < (L − 2)/4 .
Let us now explore how these two quantities approach asymptotically in the limit of M → ∞, which thus requires to take the limit of L → ∞ as well under the condition that 4M ≤ L A or M < (L − 2)/4 . To this end, here we calculate the exponents δ S (M ) and δ E (M ) by the following formulae: It is known [66] that the entanglement entropy S exact (L A ) of the exact ground state for the size L A of subsystem A is given as Equation (94) is motivated by the assumption that L A is replaced as L A ∼ M δ S , indicating that a finite M intro-duces a finite correlation length, as does a finite bond dimension in the matrix product states [64,65]. Similarly, it is also known [67][68][69] that the finite size correction to the exact ground-state energy per site is given as where ∆ε exact (L) is the energy difference between the exact ground-state energy per site for the system size L and that in the thermodynamic limit. Equation (95) is thus motivated by assuming that L ∼ M δ E . We find in Fig. 8 that these two exponents δ S (M ) and δ E (M ) approach to one in the limit of M → ∞. Thus, the asymptotic behaviors of the entanglement entropy Finally, we discuss the relation to the evolution of the single-particle orbitals in |ψ M (θ) via the local unitary time-evolution operators. As discussed in Sec. III B, the spatial extent d M of a single-particle orbital in the DQAP ansatz |ψ M (θ) with M layers of the local unitary timeevolution operators is d M = 4M + 2. Therefore, when 4M ≤ L A , for which we find that the entanglement entropy S A of |ψ M (θ) is independent of L (see Fig. 6), a single-particle orbital in |ψ M (θ) may cross one of the partitioning boundaries, but it cannot cross the other partitioning boundary. Conversely, when 4M > L A , a single-particle orbital (but not necessarily all the singleparticle orbitals) can cross both partitioning boundaries of the subsystems. Since the entanglement entropy of the free fermion system is determined by the hybridization between the two subsystems, the observation here suggests that the contribution of each partitioning boundary to the entanglement entropy S A is indeed separable in the case for 4M ≤ L A , i.e., S A ∼ S ∂AI + S ∂AII , where ∂A I(II) is the partitioning boundary and S ∂A I(II) implies the entanglement entropy from the boundary ∂A I(II) (see Apendix B).

D. Mutual information
We shall now examine the evolution of the mutual information I x,x defined in Eq. (71) for the DQAP ansatz |ψ M (θ) with increasing the number M of layers of the local unitary time-evolution operators. As discussed in Sec. II E, I x,x is a measure to quantify the entanglement between sites x and x for a quantum state. Figure 9 shows the results for the DQAP ansatz |ψ M (θ) with the variational parameters θ optimized for each M to minimize the variational energy.
We find in Fig. 9 that the mutual information I x,x is exactly zero, implying no entanglement, when |x − x | > 4M + 1. This is generally the case for any system size L. As illustrated in Fig. 10, this entanglement feature reflects the causality-cone like structure of the propagation of quantum entanglement via the local unitary timeevolution operators in the DQAC ansatz, which limits the propagation speed set by the Lieb-Robinson bound. Two causality cones for the propagation of quantum entanglement from sites x and x are indicated in Fig. 10. All the local unitary time-evolution operators inside the causality cones are connected to the origin of the cone (i.e., x or x ), while those outside the causality cones are essentially disconnected. When these two causality cones do not overlap to each other, the mutual information I x,x is zero. On the other hand, if these two causality cones overlap, we obtain I x,x = 0. We should also note that although the mutual information I x,x becomes nonzero for all values of x once the number M of layers of the local time-evolution operators satisfies L/2 ≤ 4M + 1 (see Fig. 9), more layers of the local time-evolution operators are required for the DQAP ansatz |ψ M (θ) to represent the exact ground state of the system, as discussed in Sec. III A.
The feature of the mutual information found here can also be understood on the basis of the singleparticle orbitals in the DQAP ansatz |ψ M (θ) discussed in Sec. III B. As already described in Sec. II E, the mutual information I x,x for the free fermion system is fully determined by the one-particle density matrix D {x,x } given in Eq. (72). Since ψ M (θ)|ĉ † xĉx |ψ M (θ) = 0.5 in our system, I x,x is determined solely by the off-diagonal element ψ M (θ)|ĉ † xĉx |ψ M (θ) . For example, as discussed in Sec. II E, I x,x = 0 when ψ M (θ)|ĉ † xĉx |ψ M (θ) = 0. Let us now evaluate ψ M (θ)|ĉ † xĉx |ψ M (θ) using Eq. (40). Considering that Ψ † M Ψ M = I N , we obtain that x n = 0 when |x − x | > 4M + 1, which thus also explains that I x,x = 0 when |x−x | > 4M +1 found in Fig. 9. Although I x,x becomes finite for all distances |x − x | when M satisfies 4M ≥ L/2 − 1, it is not sufficient for the DQAP ansatz |ψ M (θ) to represent the exact ground state since the single-particle orbital [Ψ M ] xn does not extend over the entire region of the system until d M = 4M + 2 ≥ L.

E. Optimized variational parameters
We shall now discuss the optimized variational parameters in the DQAP ansatz |ψ M (θ) . The natural gradient method described in Sec. II D is employed without any difficulty to optimize the variational parameters in the DQAP ansatz |ψ M (θ) , which becomes the exact ground state of the system for M = L/4 under APBC and M = (L − 2)/4 under PBC at half filling. However, we find that the optimized variational parameters are not unique and many different sets of optimized variational parameters give the same energy, as discussed in Appendix A.
Among many sets of optimum solutions for the variational parameters, we find a series of systematic solutions by gradually increasing M for a fixed system size L. Such a series is obtained as follows. We first start with a small value of M , for which the optimized variational parameters {θ can be uniquely determined for p = 1, 2. Then, we use these optimized variational parameters as the initial parameters {θ } M and optimize the variational parameters in |ψ M +1 (θ) . Here, we assumed that M is even. When M is odd, M/2 should be replaced with (M − 1)/2. With iteratively increasing M by one in this procedure, we finally obtain the series of the optimized variational parameters systematically, as shown in Fig. 11.
The characteristic features of the optimized variational parameters are summarized as follows. First of all, θ  finally decrease with increasing m when (2m − 1)/M approaches to one. This might be understood because at the last stage of the process, it would be better for the discretized quantum adiabatic process to be determined by the time-evolution operator e −iĤ f t of the final systemĤ f as in the continuous time quantum adiabatic process. To this end, the parameters θ should be small to reduce the Suzuki-Trotter decomposition error due to the discretization of time [70,71].
Notice also that there is an abrupt change of the optimized variational parameters between M = L/4 − 1 and M = L/4 (see the results for M = 39 and M = 40 in Fig. 11). As already described in Sec. III A, the optimized DQAP ansatz |ψ M (θ) with M = L/4 represents the exact ground state of the system under APBC. This abrupt change of the optimized variational parameters is associated with that of the variational energy found in Fig. 3. We also notice in Fig. 11 that the optimized variational parameters in |ψ M (θ) with M < L/4 for a given system size L converges systematically to those with M = L/4 − 1 as M increases, which are different from the optimized variational parameters in |ψ M (θ) with M = L/4. Furthermore, we find that the optimized variational parameters in |ψ M (θ) with the given number M of layers remain unchanged for different system sizes L, as long as M < L/4, which is associated with the observation that the variational energy per site, E M (L)/L, is independent of L when M < L/4, as shown in 3(b). We should also note that the optimized variational parameters in |ψ M (θ) for the system under PBC are exactly the same as those in |ψ M (θ) for the system under APBC, e.g., shown in Fig. 11, independently of the system size L, as long as M < (L − 2)/4 (also see Sec. III C). Figure 12 summarizes the optimized variational parameters in the DQAP ansatz |ψ M (θ) with M = L/4 for different system sizes L under APBC, for which |ψ M (θ) represents the exact ground state. We find that the optimized variational parameters θ = {θ } converge asymptotically to a smooth function of m for each θ (m) p (p = 1, 2) with increasing the system size L.
Finally, we examine the effective total evolution time T eff (L) of the DQAP given by where the variational parameters θ (m) p (p = 1, 2) in the DQAP ansatz |ψ M (θ) are optimized for the system size L with M = L/4 under APBC, thus representing the exact ground state, and are already shown in Fig. 12. It is highly interesting to find in Fig. 13 that the effective total evolution time T eff (L) is almost perfectly proportional to the system size L. According to the quantum adiabatic theorem, the evolution time necessary to successfully converge to the ground state of the final Hamiltonian in the continuous time quantum adiabatic process is inversely proportional to the square of the minimum energy gap during the process [35]. For the model studied here, the minimum gap appears at the final Hamiltonian, i.e., the free fermion model, and thus it is ∼ 1/L, suggesting that, according to the adiabatic theorem, the evolution time to successfully obtain the final state is proportional to L 2 . Therefore, the discretized version of the quantum adiabatic process based on the DQAP ansatz significantly improves the evolution time.

IV. IMAGINARY-TIME EVOLUTION: COMPARISON WITH THE DQAP ANSATZ
Let us now consider an ansatz inspired by the imaginary-time evolution, instead of the real-time evolution as in the DQAP ansatz |ψ M (θ) discussed in the previous sections. The imaginary-time counterpart |ϕ M (τ ) of the DQAP ansatz |ψ M (θ) defined in Eq. (27) for the free fermion system is given by where the initial state |ψ i is a product state of the local bonding states given in Eq. (25), i.e., the ground state ofV 1 , and the imaginary-time steps τ = {τ } are considered as real variational parameters that are determined so as to minimize the variational energy [72]. As in the DQAP ansatz |ψ M (θ) , |ϕ M (τ ) is constructed by repeatedly applying the local but now imaginary-time evolution operators e −τ (m) 1V 1 and e −τ (m) 2V 2 . Since these imaginary-time evolution operators are not unitary, |ϕ M (τ ) is no longer normalized.
We can follow exactly the same analysis in Sec. II for the discretized imaginary-time evolution ansatz |ϕ M (τ ) in Eq. (102) and obtain that where G M is the L × N matrix given by Here, the L × L matrices V 1 and V 2 are defined in Eqs. (83) and (84), respectively, and the L×N matrix Ψ i is given in Eq. (86). Note also that e −τ V1 and e −τ V2 are the same block diagonal matrices of e −iθV1 and e −iθV2 in Eqs. (87) and (88), respectively, except that θ is replaced with −iτ . It is now apparent that the imaginary-time evolved state |ϕ M (τ ) from a state initially prepared as a single Slater determinant state |ψ i can still be represented as a single Slater determinant state, in which each single-particle orbital is given by each column vector of G M . However, note that the imaginary-time evolved single-particle orbitals are neither normalized nor orthogonal to each other, i.e., G † M G M = I N , even though the initial single-particle orbitals are orthonormalized, i.e., The natural gradient method described in Sec. II D is straightforwardly extended to optimize the variational parameters τ = {τ } in |ϕ M (τ ) by replacing S and f in Eqs. (53) and (54), respectively, with and where F = (G † M G M ) −1 and k, k = 1, 2, · · · , 2M labeling the variational parameters {τ }. ∂ k G M is an L × N matrix defined as the first derivative of G M with respect to the kth variational parameter τ k , i.e., Figures 14(a) and 14(b) respectively show an error of the optimized variational energy and a distance defined through the fidelity as a function of the system size L for three values of M . Here, E exact (L) and |ψ exact denote the exact ground state energy and the normalized exact ground state of the system with the system size L, respectively. Although the system considered here is at the critical point where the expectation value ofĉ † xĉx decays algebraically with the distance |x − x |, we find that the variational energy of the discretized imaginary-time evolution ansatz |ϕ M (τ ) converges exponentially faster to the exact energy with increasing M . This is in sharp contrast to the results for the discretized real-time evolution ansatz |ψ M (θ) shown in Fig. 3. For example, for L < 100, we obtain the variational energy with accuracy ∆E as large as 2 × 10 −5 or less and the distance smaller than 10 −2 by using only M = 3. We should note that the exponentially fast convergence of the discretized imaginary-time evolution ansatz has also been reported for the transverse-field Ising model even at the critical point [73].
We shall now discuss how the efficiency of the discretized imaginary-time evolution ansatz |ϕ M (τ ) occurs. First of all, we should recognize that, although the discretized imaginary-time evolution ansatz |ϕ M (τ ) is composed of the local imaginary-time evolution operators, there is no limit of speed for propagating quantum entanglement via the local imaginary-time evolution operators because these local operators are non-unitary. This is indeed easily understood if we evaluate the local expectation value ϕ M (τ )|ĉ † xĉx+1 |ϕ M (τ ) . In this case, we have to treat all local imaginary-time evolution operators in |ϕ M (τ ) , no matter how far a local imaginary-time evolution operator e tτ (m) p (ĉ † yĉy+1 +ĉ † y+1ĉ y ) acting at sites y and y + 1 is distant from site x. Because of the nonunitarity, there is no cancellation of local imaginary-time evolution operators in the left and right sides of the local expectation value ϕ M (τ )|ĉ † xĉx+1 |ϕ M (τ ) , implying that there is no causality-cone like structure for the propagation of quantum entanglement illustrated in Fig. 4.
Second, for the free fermion system, we can understand the efficiency of the discretized imaginary-time evolution ansatz |ϕ M (τ ) in terms of the imaginary-time evolution of the single-particle orbitals in G M . For this purpose, we introduce the following L × N matrix: for m = 0, 1, · · · , M with g 0 = Ψ i and g M = G M to represent the single-particle orbitals at an intermediate imaginary time. Here, the variational parameters τ = {τ } are optimized for |ϕ M (τ ) to minimize the variational energy and these parameters are used to define the matrix g m in Eq. (110). Note that g m is an analog to Φ m introduced in Eq. (89) for the DQAP ansatz |ψ M (θ) . Figures 15(a)-15(c) show the discretized imaginarytime evolution of all matrix elements in three matrices g 0 , g 1 , and g 2 = G 2 when M = 2. As in the case of a single-particle orbital in Φ m , we can readily find that a single-particle orbital in g m extends spatially four lattice space every time the imaginary time m increases by one: the spatial extentd m (in unit of lattice constant) of a single-particle orbital in g m isd m = 4m + 2, exactly the same as the spatial extent d m of a single-particle orbital in Φ m for the real-time evolution (see Sec. III B). This is simply because the imaginary-and real-time evolutions are both governed by the spatially local evolution operators. Therefore, as in the case of the real-time evolution ansatz |ψ M (θ) , the single-particle orbitals in g m can extend over the entire system only whend m reaches to the system size L.
However, unlike the case of the real-time evolution ansatz, the single-particle orbitals in g m are not orthonormalized, i.e., g † m g m = I N for m > 0, as shown in Fig. 15(d). As a result, (g † m g m ) −1 becomes non-local in the sense that [(g † m g m ) −1 ] xx = 0 even when sites x and x are distant from each other [see Fig. 15(e)], although g † m g m might be local. This has a significant consequence when we evaluate the expectation value ϕ M (τ )|ĉ † xĉx |ϕ M (τ ) , taking also into account the normalization of |ϕ M (τ ) . Using Eq. (40), we can show that is also non-local even for M L, as shown in Fig. 15(f), implying that the expectation value ofĉ † xĉx is non-zero even when sites x and x are far apart. This should be contrasted with the case of the the DQAP ansatz |ψ M (θ) , where the corresponding expectation value is given as (112) because the real-time evolution ansatz |ψ M (θ) is normalized, i.e., Ψ † M Ψ M = I N , and thus the expectation value ofĉ † xĉx is zero when sites x and x are far apart, provided that M is not large enough, as discussed in Sec. III D.
Consequently, the discretized imaginary-time evolution ansatz |ϕ M (τ ) acquires the global correlation with the extremely small number of M . Figure 16 shows the (a) g 0 mutual information I x,x of the optimized |ϕ M (τ ) with M = 1, 2, 3. We find that, although I x,x for M = 1 shows exponential decay as a function of distance |x−x |, it is drastically improved with increasing M and the mutual information I x,x evaluated for M = 2 already almost coincides with the exact value, despite that there are only four variational parameters for M = 2. We should also emphasize that no causality-cone like structure is observed in the mutual information I x,x of the discretized imaginary-time evolution ansatz |ϕ M (τ ) , which is in sharp contrast to the results for the discretized real-time evolution ansatz |ψ M (θ) shown in Fig. 9.
where the variational parameters {τ (m) p } are optimized for each M to minimize the variational energy. First, it is noticed thatβ(L) exhibits the system size dependence, which is different from that found for the discretized realtime evolution ansatz |ψ M (θ) shown in Fig. 13. It is also important to notice thatβ(L) is not proportional to the system size L. On one hand, one would expect that a large τ (m) p is preferable to reach the ground state fast, i.e., with less number of variational parameters, in a standard sense of the imaginary-time evolution. On the other hand, a large τ (m) p might introduce bias in approximating the continuous imaginary-time evolution by the discretized evolution via the Suzuki-Trotter decomposition [70,71]. Therefore, the optimized solution should be determined by compromising these two factors. The discretized imaginary-time evolution ansatz |ψ M (θ) finds the best solution available within a given value of M .

V. SUMMARY AND DISCUSSION
As a quantum-classical hybrid algorithm to generate a desired quantum state in a quantum circuit, we have studied the DQAP ansatz |ψ M (θ) to represent the ground state of the one-dimensional free fermion system. The DQAP ansatz |ψ M (θ) considered here is inspired by the QAOA and is composed of M layers of two elementary sets of local time-evolution operators acting on neighboring sites (i.e., quabit), as illustrated in Fig. 1. By numerically optimizing the variational parameters } so as to minimize the variational energy, we have found that the exact ground state can be attained by the DQAP ansatz |ψ M (θ) with the number M B of layers as large as (L − 2)/4 for PBC and L/4 for APBC, i.e., the minimum number of M set by the Lieb-Robinson bound for the propagation of quantum entanglement via the local time-evolution operators (see Fig. 4). Our results thus suggest that the DQAP ansatz |ψ M (θ) is the ideal ansatz to represent the exact ground state based on the quantum adiabatic process. Indeed, in the DQAP scheme, the exact ground state is prepared by a shallowest quantum circuit with linear depth, containing O(L 2 ) single-qubit and CNOT gates, where L is the number of sites in the system, i.e., the number of qubits.
We have also found that the optimized DQAP ansatz |ψ M (θ) with M less than M B exhibits another series of the states that are independent of the system size L. We have shown that the entanglement entropy S A of subsystem A and the variational energy E M (L)/L per site evaluated for these states with 4M ≤ L A and M < M B , respectively, fall into smooth universal functions of M , independently of the system size L and the boundary conditions. This implies that the entanglement acquired by the DQAP ansatz |ψ M (θ) with a finite M is bounded, as in the case of the matrix product states with a finite bond dimension [39,40]. Moreover, we have found that the entanglement entropy S A and the energy difference between the variational energy and the exact one ∆ε = E M (L)/L − ε ∞ behave asymptotically as S A ≈ 1 3 ln M and ∆ε ∼ M −2 , respectively [64,65]. We have also analyzed the evolution of the singleparticle orbitals in the DQAP ansatz |ψ M (θ) via the local time-evolution operators and the mutual information of |ψ M (θ) to explore how quantum entanglement is developed in the quantum state. The latter quantity also reveals the causality-cone like structure of the propagation of quantum entanglement via the local time-evolution operators. Furthermore, we have found that the optimized variational parameters θ = {θ } in the DQAP ansatz |ψ M (θ) converge to a smooth function of m for each θ (m) p (p = 1, 2), which is quite different from the linear scheduling functions expected when the quantum adiabatic process is naively discretized in time. We have also found that the effective total evolution time T eff (L) of the optimized variational parameters θ in the DQAP ansatz |ψ M (θ) with M = M B , thus representing the exact ground state, scales linearly with the system size L, as opposed to L 2 expected in the continuous-time quantum adiabatic process.
We have also explored the discretized imaginary-time evolution ansatz |ϕ M (τ ) , an imaginary-time counterpart of the DQAP ansatz |ψ M (θ) , for the same free fermion system. Similarly to the DQAP ansatz |ψ M (θ) , the discretized imaginary-time evolution ansatz |ϕ M (τ ) is composed of M layers of two elementary sets of local imaginary-time evolution operators acting on neighboring sites. We have found that the convergence of |ϕ M (τ ) to the exact ground state is exponentially fast in terms of the number M of layers, as compared to that of the DQAP ansatz |ψ M (θ) , although both ansatze are composed of the local evolution operators. This difference is attributed to the fact that the imaginary-time evolution operator is not unitary and thus there is no limit of speed for the propagation of quantum entanglement via the local non-unitary imaginary-time evolution operators. In particular, for the free fermion system, we can show directly that the expectation value ϕ M (τ )|ĉ † xĉx |ϕ M (τ ) / ϕ M (τ )|ϕ M (τ ) is non-local even when M L, regardless of the distance |x − x |, because the discretized imaginary-time evolution ansatz |ϕ M (τ ) is not normalized, i.e., ϕ M (τ )|ϕ M (τ ) = 1. This is in sharp contrast to the case of the DQAP ansatz |ψ M (θ) , in which the corresponding expectation value ψ M (θ)|ĉ † xĉx |ψ M (θ) is zero when the two causality cones set by the Lieb-Robinson bound formed from the origins at sites x and x do not overlap (see Fig. 10). Our result thus implies that if the local non-unitary imaginary-time evolution operator can be implemented in a quantum circuit by using O(1) local single-and twoqubit unitary gates, one can prepare the ground state in this scheme by a much shallower quantum circuit with depth O(1). However, it is challenging to represent a local non-unitary operator by O(1) local single-and twoqubit unitary gates, specially, for a quantum state at criticality. The free fermion system considered here is at the critical point where the correlation function ĉ † xĉx decays algebraically with |x − x | and thus the correlation is extended over the entire system. In a critical system, it is intuitively understood that at least L 2 local unitary twoqubit gates are required in a quantum circuit to represent the quantum entanglement of the state for the system with L sites. Therefore, also in this sense, the DQAP ansatz |ψ M (θ) is an ideally compact form to represent the ground state of this system. However, it is not trivial for more general cases such as an interacting fermion system. It is thus valuable to consider a possible improvement in the quantum adiabatic process, for example, by introducing navigation proposed in the VanQver algorithm [74], for reducing the complexity of quantum processes. It is also an interesting extension to introduce a non-unitary process by inserting measurements during the quantum adiabatic process [75,76].
We have also found that the natural gradient method can optimize the variational parameters in the DQAP ansatz |ψ M (θ) without any difficulty. Even with randomly chosen initial variational parameters, the opti-mization method can eventually find sets of optimized variational parameters to converge to the lowest variational energy (see Appendix A), implying that there is no problem such as the barren plateaus phenomena [7]. However, this could be due to the fact that for the free fermion system, the independent matrix elements in the DQAP ansatz |ψ M (θ) can be significantly reduced (see discussion in Sec. III B). Therefore, it is desirable to examine the case for an interacting fermion system, for example. This is left for a future study.
The focus in this paper was limited on the free fermion system, where a time-evolved N -fermion state can still be described by a single Slater determinant state, and therefore any quantum advantage is expected in simulating this system on a quantum computer. However, this system is one of the ideal systems to test the operations of NISQ devices because the quantum state described by the DQAP ansatz |ψ M (θ) is highly entangled but can still be treated in large systems on a classical computer.
are clearly different. Consequently, the single-particle orbitals in the optimized DQAP ansatz |ψ M (θ) have different shapes, depending on the variational parameters θ (m) p . However, this does not alter the conclusion discussed in Sec. III B, qualitatively.
Appendix B: Boundary contribution to entanglement entropy As described in Sec. III C, the entanglement entropy S A of the optimized DQAP ansatz |ψ M (θ) is determined by the number M of layers in the DQAP ansatz, independently of the system size L and the boundary conditions, provided that 4M ≤ L A and L A ≤ LĀ, where L A (LĀ) is the size of subsystem A (complement of subsystem A) and L = L A + LĀ. This finding suggests that the entanglement entropy S A is separable to the contributions from the partitioning boundaries ∂A I and ∂A II between the two subsystems, i.e., S A ∼ S ∂AI + S ∂AII , where S ∂A I(II) implies the entanglement entropy from the boundary ∂A I(II) . Note that the partitioning boundaries are assumed not to break any local bonding state in the initial state |ψ i (see Fig. 7). Here, we discuss more details of this point through the one-particle density matrix D A defined in Eq. (60).
Let us first introduce the one-particle density matrix D U of the whole system U = A +Ā as whereĉ * (ĉ t ) is the matrix transpose ofĉ † (ĉ) in Eq. (3) and Ψ M is given in Eq. (81). We assume that A = {1, 2, 3, · · · , L A } andĀ = {L A +1, L A +2, L A +3, · · · , L}, for simplicity. Then, we can readily show that because Ψ † M Ψ M = I N . This implies that the eigenvalues of D U are either 0 or 1.
Let us now write D U as where D AA is an L A × L A matrix, corresponding to the one-particle density matrix D A of the subsystem A defined in Eq. (60), and DĀĀ is an LĀ × LĀ matrix. Due to the idempotence of D U in Eq. (B3), we find that irrespectively of the values of the variational parameters in the DQAP ansatz |ψ M (θ) . Here, the spatial extent d M = 4M + 2 of the single-particle orbitals in Ψ M suggests that D U = Ψ * M Ψ t M is a band-like matrix with 8M + 2 non-zero elements in each row and each column, while D AĀ has an opposite matrix structure with nonzero elements appearing at the upper right and lower left corners, and the apparent rank of D AĀ is 8M . For example, when L = 16, N = 8, and M = 1, the matrix structures of Ψ * M Ψ t M = D U are schematically given as   * * * * * * 0 0 0 0 0 0 * * * * * * * * * * 0 0 0 0 0 0 * * * * * * * * * * * * 0 0 0 0 0 0 * * * * * * * * * * 0 0 0 0 0 0 * * * * * * * * * * * * 0 0 0 0 0 0 * * * * * * * * * * 0 0 0 0 0 0 0 0 * * * * * * * * * * 0 0 0 0 0 0 * * * * * * * * * * 0 0 0 0 0 0 0 0 * * * * * * * * * * 0 0 0 0 0 0 * * * * * * * * * * 0 0 0 0 0 0 0 0 * * * * * * * * * * 0 0 0 0 0 0 * * * * * * * * * * * * 0 0 0 0 0 0 * * * * * * * * * * 0 0 0 0 0 0 * * * * * * * * * * * * 0 0 0 0 0 0 * * * * * * * * * * 0 0 0 0 0 0 * * * * * * where " * " indicates a non-zero element and D AĀ is the upper right quadrant of the matrix in the right hand side. However, due to the characteristic structure of Ψ M , a single-particle orbital being extended spatially by two lattice space in each spatial direction every time the local time-evolution operators are applied, we find that the Gaussian elimination eliminates a half of the non-zero  Fig. 7(a), there are L A /2 − 2M singleparticle orbitals in the DQAP ansatz |ψ M (θ) that do not cross either side of the partitioning boundaries between the two subsystems and stay inside the subsystem A. These L A /2 − 2M single-particle orbitals contribute to the eigenvalues of D A with δ l = 1. Due to the particlehole symmetry, the eigenvalues of D A should appear symmetrically around 1/2. Therefore, there exist L A /2−2M eigenvalues of D A with δ l = 0, corresponding to the unoccupied single-particle orbitals that stay inside the subsystem A without crossing the partitioning boundaries. The remaining 4M eigenvalues of D A are neither 0 nor 1, i.e., 0 < δ l < 1. These contributions are due to the singleparticle orbitals (2M ) in the DQAP ansatz |ψ M (θ) that cross either side of the partitioning boundaries between the two subsystems and the hole counterparts (2M ) due to the particle-hole symmetry. These eigenvalues of D A with 0 < δ l < 1 contribute to the non-zero entanglement entropy S A in Eq. (69). Figure 18 shows  Fig. 18 are also plotted. M = 5 and two different system sizes L = 40 and 80, assuming that L A = L/2 and thus 4M ≤ L A . First, we can notice that the eigenvalues are all symmetric around 1/2, as expected due to the particle-hole symmetry. Second, we can confirm that all the eigenvalues are neither 0 nor 1 for L = 40. Third, there are L A /2 − 2M = 10 eigenvalues with δ l = 1 as well as δ l = 0 for L = 80. In addition, we can find that other eigenvalues different from 0 or 1 for L = 80 are identical to the eigenvalues of D A found for L = 40. These eigenvalues δ l with 0 < δ l < 1 contribute to the entanglement entropy S A , and therefore this finding is in good agreement with the result that the entanglement entropy S A of the optimized DQAP ansatz |ψ M (θ) with 4M ≤ L A is independent of the system size L (see Fig. 6). Fourth, the eigenvalues δ l with 0 < δ l < 1 are pairwise degenerate. Considering that only these 4M eigenvalues δ l with 0 < δ l < 1 contribute to the entanglement entropy S A and these eigenvalues correspond to the single-particle orbitals in the DQAP ansatz |ψ M (θ) that cross either side of the partitioning boundaries and the hole counterparts, we interpret the pairwise degeneracy as a fingerprint that the contribution to the entanglement entropy S A can be separated from the two partitioning boundaries ∂A I and ∂A II , i.e, S A ∼ S ∂AI + S ∂AII . Indeed, as shown in Fig. 19, the pairwise degeneracy disappears once 4M > L A [77].