update en doc for algorithm

docs/sphinx/en/algorithm/algorithms.rst

@@ -170,6 +170,103 @@ The singular values obtained from the SVD of the matrix are used as the mean fie
 
 Although the computation cost of the simple update is cheaper than that of the full update, it is known that the simple update shows strong initial state dependence and it tends to overestimate the local magnetization. Thus, for complicated problems, we need to carefully check results obtained by the simple update. 
 
+Real-time evolution by iTPS
+===========================
+The algorithms of imaginary time evolution used for computing the ground state, such as the simple update method and the full update method, can also be used to calculate the real-time evolution of a quantum state. In TeNeS, similarly to the case of imaginary time evolution, the quantum state at time :math:`t`
+
+.. math::
+   |\Psi(t)\rangle = e^{-it\mathcal{H}} |\Psi_0\rangle,
+
+is approximated by iTPS, which allows for the calculation of approximate time evolution. The difference between imaginary and real-time evolution lies only in whether the coefficient of the Hamiltonian :math:`\mathcal{H}` in the exponent is :math:`-\tau` or :math:`-it`, hence real-time evolution can also be computed using the same simple update and full update methods applied in imaginary time evolution, by employing the Suzuki-Trotter decomposition.
+
+Real-time evolution using iTPS (and other tensor network states) differs significantly from imaginary time evolution used for ground state calculation in two main aspects.
+
+One major difference is the size of the quantum entanglement of the target quantum state. In imaginary time evolution, as the evolution progresses towards the ground state, the quantum entanglement of the state does not become excessively large. Thus, the description by iTPS works well. However, in real-time evolution, typically (unless the initial state's iTPS is an eigenstate of the Hamiltonian), quantum entanglement can increase over time. To maintain the approximation accuracy of iTPS, it is necessary to increase the bond dimension of iTPS as the time gets longer. Naturally, increasing the bond dimension also increases computational costs, so with realistic computational resources, accurately approximating real-time evolution using iTPS is limited to short times. The applicable time range depends on the model, but for example, in spin models, the limit is often around a time :math:`t = O(1/J)` with respect to the typical interaction strength :math:`J`.
+
+Another difference is the characteristics of the physical phenomenon to be reproduced. When using imaginary time evolution to calculate the ground state, it is sufficient to reach the ground state after a sufficiently long evolution, so minor deviations from the correct path of imaginary time evolution are not a significant issue. On the other hand, in real-time evolution, there is often interest not only in the final state but also in the time evolution of the quantum state itself. To accurately approximate the path of time evolution, it is necessary to not only increase the bond dimension of iTPS but also to make the time increment :math:`\delta t` of the Suzuki-Trotter decomposition sufficiently small. Depending on the situation, it may be more efficient to use higher-order Suzuki-Trotter decompositions. In TeNeS, it is possible to handle higher-order Suzuki-Trotter decompositions by editing the ``evolution`` section of the input file that is ultimately entered into tenes.
+
+Finite temperature simulation
+===========================
+So far, we considered the tensor network representation of a pure state :math:`|\Psi\rangle`, but similarly, we can consider the tensor network representation for a mixed state at finite temperature
+
+.. math::
+   \rho(\beta) = \frac{e^{-\beta \mathcal{H}}}{\mathrm{Tr} e^{-\beta \mathcal{H}}}
+
+where :math:`\beta` represents the inverse temperature corresponding to temperature :math:`T` as :math:`\beta = 1/T`.
+
+Similarly to pure states, if we consider a system of :math:`N` quantum spins with :math:`S=1/2` at finite temperature, the mixed state can be expressed as
+
+.. math::
+   \rho(\beta) = \sum_{s_i=\uparrow, \downarrow, s_i' = \uparrow, \downarrow} \left(\rho(\beta)\right)_{s_1,s_2,\dots, s_N}^{s_1', s_2', \dots, s_N'} |s_1', s_2', \dots, s_N'\rangle \langle s_1, s_2, \dots, s_N|
+
+The expansion coefficients :math:`\left(\rho(\beta)\right)_{s_1,s_2,\dots, s_N}^{s_1', s_2', \dots, s_N'}` can be expressed, for example, using a Matrix Product Operator (MPO), generalized from MPS to matrices (operators), as
+
+.. math::
+   \left(\rho^{\mathrm{MPO}}(\beta)\right)_{s_1,s_2,\dots, s_N}^{s_1', s_2', \dots, s_N'}  = T^{(1)}[s_1, s_1']T^{(2)}[s_2, s_2']\cdots T^{(N)}[s_N, s_N']
+
+and the corresponding diagram can be drawn as
+
+.. image:: ../../img/MPO.*
+   :align: center
+
+For mixed states with translational symmetry, just like in the case of pure states, an infinite MPO (iMPO) can represent the state of an infinite system by repeating the same tensor infinitely. For example, for a one-dimensional, two-site translational symmetric state, the corresponding iMPO diagram would be
+
+.. image:: ../../img/iMPO.*
+   :align: center
+
+As a tensor network represention of mixed states, in TeNeS, we handle a two-dimensional infinite tensor product operator (iIPO) :ref:`[TPO] <Ref-TPO>`, specifically assuming a square lattice network with translational symmetry. The diagram for such an iTPO can be written as
+
+.. image:: ../../img/iTPO.*
+   :align: center
+
+In TeNeS, the mixed state at finite temperature :math:`\rho(\beta)` is computed using imaginary time evolution from the initial state corresponding to infinite temperature :math:`\rho(\beta=0)`
+
+.. math::
+    \rho(\beta) = e^{-\frac{\beta}{2} \mathcal{H}} \rho(0) e^{-\frac{\beta}{2} \mathcal{H}}
+
+Note that at infinite temperature, the density matrix is the identity matrix. From this property, for example, the iMPO representation of the state at infinite temperature becomes a tensor product of local identity matrices, and the diagram in this case would be drawn as
+
+.. image:: ../../img/iMPO_T0.*
+   :align: center
+
+with "lines" corresponding to the local identity matrix.
+
+The imaginary time evolution of a mixed state is calculated by a simple extension of the imaginary time evolution for pure states, as an approximate imaginary time evolution within the iTPO representation. The Suzuki-Trotter decomposition, simple update method, and full update method used for pure states can be almost directly applied to the case of mixed states. (TeNeS does not support the full update currently.)
+
+The local minimization problem for mixed states can be described as
+
+.. math::
+   \min \left \Vert \rho_{\tau}^{\mathrm{iTPO}} - e^{-\frac{\tau}{2} \mathcal{H}_{ij}/2} \rho^{\mathrm{iTPO}} e^{-\frac{\tau}{2} \mathcal{H}_{ij}}\right \Vert^2
+
+and the corresponding diagram, for clarity in the form of iMPO, would be 
+
+.. image:: ../../img/iMPO_ITE_local.*
+   :align: center
+
+The biggest difference between the computations of finite temperature states by iTPO and pure states by iTPS appears in the tensor network for expectation value calculations. The expectation value of a physical quantity :math:`O` for a given mixed state :math:`\rho` is calculated as
+
+.. math::
+   \langle O \rangle_\rho = \frac{\mathrm{Tr} (\rho O)}{\mathrm{Tr} \rho}.
+
+The trace :math:`\mathrm{Tr}` corresponds to connecting the corresponding upper and lower legs of the iTPO. Using a tensor obtained by connecting upper and lower legs of a local tensor in iTPO,
+
+.. image:: ../../img/trace_tensor.*
+   :align: center
+
+the denominator :math:`\mathrm{Tr} \rho` becomes the same structure as the two-dimensional square lattice diagram appeared in the expectation values for pure states. Thus, we can apply the same approximate calculation using corner transfer matrix representation and CTMRG.
+
+The computation cost of CTMRG for the corner transfer matrix representation with bond dimension :math:`\chi` and iTPO with bond dimension :math:`D` scales with :math:`O(\chi^2 D^4)` and :math:`O(\chi^3 D^3)`. Note that this computation cost is smaller compared to CTMRG for pure states with the same bond dimension :math:`D`. The difference is due to the bond dimension of the tensor indicated by the black circle being :math:`D^2` in pure state calculations, while :math:`D` for mixed states. Correspondingly, the bond dimension :math:`\chi` of the corner transfer matrices can be increased proportionally to :math:`D`, i.e., :math:`\chi \propto O(D)`. Under this condition, the computation cost of CTMRG becomes :math:`O(D^6)`, and the required memory amount becomes :math:`O(D^4)`. Thus, the computation cost of finite temperature calculations using iTPO is significantly lower than that of iTPS with the same :math:`D`. It allows us to use larger bond dimensions :math:`D` in finite temperature calculations.
+
+Similarly to pure states, once the converged corner transfer matrices and edge tensors are computed, :math:`\mathrm{Tr} (\rho O)` can also be efficiently calculated. For example, when we define the tensor containing the operator as
+
+.. image:: ../../img/trace_Sz.*
+   :align: center
+
+, the local magnetization :math:`\mathrm{Tr} (\rho S_i^z)` is calculated using the same diagram as :math:`\langle \Psi|S_i^z|\Psi\rangle`.
+
+Lastly, it is important to mention the drawbacks of approximation by iTPO. The density matrix of a mixed state is Hermitian and positive semidefinite, with non-negative eigenvalues. However, when approximating the density matrix with iTPO, this positive semidefiniteness is not guaranteed, and physical quantities calculated from the iTPO approximation might exhibit unphysical behavior, such as energies lower than the ground state energy. This is a problem of iTPO representation, and cannot be avoided just by improving the accuracy of CTMRG in expectation value calculation by increasing the bond dimension :math:`\chi`. To recover physical behavior, it is necessary to increase the bond dimension :math:`D` of iTPO to improve the approximation accuracy of the density matrix.
+
+As an alternative representation to avoid such unphysical behavior, a method has been proposed using purification of the density matrix, representing the purified density matrix with iTPO :ref:`[Purification] <Ref-Purification>`. However, in this case, the diagram appearing in the expectation value calculation becomes a double-layer structure similar to pure states. This structre requires a larger computational cost, and the manageable bond dimension :math:`D` becomes smaller than in the direct iTPO representation.
 
 
 .. rubric:: References
@@ -203,3 +300,12 @@ H. G. Jiang *et al.*, *Accurate Determination of Tensor Network State of Quantum
 
 [QR]
 L. Wang *et al.*, *Monte Carlo simulation with tensor network states*, Phys. Rev. B **83**, 134421 (2011). `link <https://doi.org/10.1103/PhysRevB.83.134421>`__
+
+.. _Ref-TPO:
+
+[TPO]
+A. Kshetrimayum, M. Rizzi, J. Eisert, and R. Orús, *Tensor Network Annealing Algorithm for Two-Dimensional Thermal States*, Phys. Rev. Lett. **122**, 070502 (2019). `link <https://doi.org/10.1103/PhysRevLett.122.070502>`__
+
+.. _Ref-Purification:
+[Purification]
+P. Czarnik, J. Dziarmaga, and P. Corboz, *Time evolution of an infinite projected entangled pair state: An efficient algorithm*, Phys. Rev. B **99**, 035115 (2019). `link <https://doi.org/10.1103/PhysRevB.99.035115>`__; P. Czarnik and J. Dziarmaga, *Time evolution of an infinite projected entangled pair state: An algorithm from first principles*, Phys. Rev. B **98**, 045110 (2018). `link <https://doi.org/10.1103/PhysRevB.98.045110>`__