Math fixes

README.md

@@ -5,7 +5,7 @@
 The code uses matrix-product-state (MPS) and matrix-product-operator (MPO) data structures, implemented in C++ (using the [ITensor](https://itensor.org/) library) for maximizing the performance with multithreaded computations, and wrapped by a Python layer with an easy-to-use interface and rich plotting features. The solver integrates in fixed time steps the Markovian master equation for the density matrix,
 
   ```math
-\frac{\partial}{\partial t}\rho = -\frac{i}{\hbar}[H,\rho]+\mathcal{D}[\rho]
+\frac{\partial}{\partial t}\rho = -\frac{i}{\hbar}[H,\rho]+\mathcal{D}[\rho].
   ```
 
 The solver supports various initial states and observables, with the Hamiltonian and dissipator terms having time-independent coefficients;

docs/API_DOCS.md

@@ -17,61 +17,56 @@ All the model parameters are passed through the parameters dictionary in the cla
 The default value for each parameter is indicated with an equality sign (only the first three parameters below have no default value and must be specified when the solver is invoked using the Python interface).
 
 * Basic parameters:
-    * N (int): The number of qubits in the lattice. This solver requires <img src="https://render.githubusercontent.com/render/math?math=N \geq 2" style="vertical-align:bottom">.
-    * t_final (float): The final simulation time (<img src="https://render.githubusercontent.com/render/math?math=t_f" align=middle>).
-    * tau (float): The discrete time step <img src="https://render.githubusercontent.com/render/math?math=\tau" style="vertical-align:bottom"> used in the time evolution.
-    * t_init = 0 (float): The initial simulation time (<img src="https://render.githubusercontent.com/render/math?math=t_0" align=middle>). Must obey <img src="https://render.githubusercontent.com/render/math?math=t_0 \le t_f" align=middle>.
+    * N (int): The number of qubits in the lattice. This solver requires $N \geq 2$.
+    * t_final (float): The final simulation time ($t_f$).
+    * tau (float): The discrete time step $\tau$ used in the time evolution.
+    * t_init = 0 (float): The initial simulation time ($t_0$). Must obey $t_0 \le t_f$.
     * output_files_prefix = "lindblad" (str): The path and file name prefix to be used for the input file generated for the solver, as well as output files generated by the solver. The default value indicates saving files in the working directory with the file name prefix "lindblad". If only a directory path is given, the "lindblad" file name prefix will be added.
     * b_unique_id = False (bool): If True, a unique id will be generated for the simulation and appended to all generated file name prefixes.
 * Hamiltonian coefficients:
-    * h_x = 0 (float): The <img src="https://render.githubusercontent.com/render/math?math=h_{x,i}" > coefficient in the Hamiltonian. If a vector is given, it specifies <img src="https://render.githubusercontent.com/render/math?math=h_{x,i}" style="vertical-align:middle"> for each qubit. If a scalar is given, it is uniform for all qubits.
-    * h_y = 0 (float): The <img src="https://render.githubusercontent.com/render/math?math=h_{y,i}" style="vertical-align:bottom"> coefficient in the Hamiltonian. If a vector is given, it specifies <img src="https://render.githubusercontent.com/render/math?math=h_{y,i}" style="vertical-align:bottom"> for each qubit. If a scalar is given, it is uniform for all qubits.
-    * h_z = 0 (float): The <img src="https://render.githubusercontent.com/render/math?math=h_{z,i}" style="vertical-align:bottom"> coefficient in the Hamiltonian. If a vector is given, it specifies <img src="https://render.githubusercontent.com/render/math?math=h_{z,i}" style="vertical-align:bottom"> for each qubit. If a scalar is given, it is uniform for all qubits.
-    * J_z = 0 (float): The <img src="https://render.githubusercontent.com/render/math?math=J^z_{ij}" style="vertical-align:bottom"> coefficient in the interaction part of the Hamiltonian. If a matrix is given, it specifies <img src="https://render.githubusercontent.com/render/math?math=J^z_{ij}" style="vertical-align:bottom"> for each pair of qubits. If a scalar is given, it is uniform for all qubits of a lattice. If either one of <img src="https://render.githubusercontent.com/render/math?math=J" style="vertical-align:bottom"> or <img src="https://render.githubusercontent.com/render/math?math=J_z" style="vertical-align:bottom"> is a matrix, then the other one must be either a matrix as well, or 0.
-    * J = 0 (float): The <img src="https://render.githubusercontent.com/render/math?math=J_{ij}" style="vertical-align:bottom"> coefficient in the interaction part of the Hamiltonian. If a matrix is given, it specifies <img src="https://render.githubusercontent.com/render/math?math=J_{ij}" style="vertical-align:bottom"> for each pair of qubits. If a scalar is given, it is uniform for all qubits of a lattice. If either one of <img src="https://render.githubusercontent.com/render/math?math=J" style="vertical-align:bottom"> or <img src="https://render.githubusercontent.com/render/math?math=J_z" style="vertical-align:bottom"> is a matrix, then the other one must be either a matrix as well, or 0.
+    * h_x = 0 (float): The $h_{x,i}$ coefficient in the Hamiltonian. If a vector is given, it specifies $h_{x,i}$ for each qubit. If a scalar is given, it is uniform for all qubits.
+    * h_y = 0 (float): The $h_{y,i}$ coefficient in the Hamiltonian. If a vector is given, it specifies $h_{y,i}$ for each qubit. If a scalar is given, it is uniform for all qubits.
+    * h_z = 0 (float): The $h_{z,i}$ coefficient in the Hamiltonian. If a vector is given, it specifies $h_{z,i}$ for each qubit. If a scalar is given, it is uniform for all qubits.
+    * J_z = 0 (float): The $J^z_{ij}$ coefficient in the interaction part of the Hamiltonian. If a matrix is given, it specifies $J^z_{ij}$ for each pair of qubits. If a scalar is given, it is uniform for all qubits of a lattice. If either one of $J$ or $J_z$ is a matrix, then the other one must be either a matrix as well, or 0.
+    * J = 0 (float): The $J_{ij}$ coefficient in the interaction part of the Hamiltonian. If a matrix is given, it specifies $J_{ij}$ for each pair of qubits. If a scalar is given, it is uniform for all qubits of a lattice. If either one of $J$ or $J_z$ is a matrix, then the other one must be either a matrix as well, or 0.
 * Dissipation coefficients:
-    * g_0 = 0 (float): The <img src="https://render.githubusercontent.com/render/math?math=g_{0,i}" style="vertical-align:bottom"> coefficient in the Lindbladian. If a vector is given, it specifies <img src="https://render.githubusercontent.com/render/math?math=g_{0,i}" style="vertical-align:bottom"> for each qubit. If a scalar is given, it is uniform for all qubits.
-    * g_1 = 0 (float): The <img src="https://render.githubusercontent.com/render/math?math=g_{1,i}" style="vertical-align:bottom"> coefficient in the Lindbladian. If a vector is given, it specifies <img src="https://render.githubusercontent.com/render/math?math=g_{1,i}" style="vertical-align:bottom"> for each qubit. If a scalar is given, it is uniform for all qubits.
-    * g_2 = 0 (float): The <img src="https://render.githubusercontent.com/render/math?math=g_{2,i}" style="vertical-align:bottom"> coefficient in the Lindbladian. If a vector is given, it specifies <img src="https://render.githubusercontent.com/render/math?math=g_{2,i}" style="vertical-align:bottom"> for each qubit. If a scalar is given, it is uniform for all qubits.
+    * g_0 = 0 (float): $g_{0,i}$ coefficient in the Lindbladian. If a vector is given, it specifies $g_{0,i}$ for each qubit. If a scalar is given, it is uniform for all qubits.
+    * g_1 = 0 (float): The $g_{1,i}$ coefficient in the Lindbladian. If a vector is given, it specifies $g_{1,i}$ for each qubit. If a scalar is given, it is uniform for all qubits.
+    * g_2 = 0 (float): The $g_{2,i}$ coefficient in the Lindbladian. If a vector is given, it specifies $g_{2,i}$ for each qubit. If a scalar is given, it is uniform for all qubits.
 * Initial state:
     The initial state of the simulator either takes a predefined default value, is defined explicitly by the parameters, or is loaded from files saved by a previous solver execution. Using the initialization parameters the state can be initialized to a product state as detailed below, and then optional controlled-Z (CZ) gates can be applied to (some) pairs of qubits, which allows creating an initial (possibly nonideal) graph state.
 
     * init_product_state = "+z": A float (or a one-element float tuple), a tuple of two or three floats, a two-character string, or a length-N vector of such entries. A single entry applies to all qubits.
-         * The supported strings take either the form "<img src="https://render.githubusercontent.com/render/math?math=\pm a" style="vertical-align:bottom">", or "id". The latter indicates the fully mixed-state (identity) density matrix, and the former indicate an eigenstate of one of the Pauli matrices,
-
-           <img src="https://render.githubusercontent.com/render/math?math=\sigma_i^a\left|\pm a_i\rangle = \pm \right|\pm a_i\rangle." style="vertical-align:bottom">
-
+         * The supported strings take either the form "+a", "-a", or "id". The latter indicates the fully mixed-state (identity) density matrix, and the former indicate an eigenstate of one of the Pauli matrices, $\sigma_i^a\left|\pm a_i\rangle = \pm \right|\pm a_i\rangle.$
          * If an entry is a float, it should be between 0 and 1, and indicates the population of the "|0>" ("|up>" state, aligned with +z) in a diagonal mixed-state density matrix for the qubit.
          * If an entry is a tuple of two floats (theta, phi), it indicates an arbitrary pure-state superposition in Bloch-sphere polar coordinates with the formula "cos(theta/2)|0> + sin(theta/2) exp{i phi}|1>".
          * If an entry is a tuple of three floats (a, b, c), it indicates the density matrix "a|0><0| + (1-a)|1><1| + [(b + ic) |0><1| + H.c.]"
 
 	* init_cz_gates = []. A list of integer tuples that specify the qubit pairs for performing a controlled-Z gate on, after they are initialized according to 'init_product_state'.
-    * init_graph_state = []. A list of integer tuples that specify the qubit pairs for performing a controlled-Z gate on, to generate an initial graph state, after first initializing all along the +x axis. If 'init_graph_state' is used, all other initialization parameters should be left empty. The qubit pairs will represent j and k in the graph state formula;
-
-      <img src="https://render.githubusercontent.com/render/math?math=\left|\psi_0\rangle=%20\prod_{(j,k)\in%20V}{CZ}[j,k]%20\prod_i%20\right|%2b%20x_i\rangle" style="vertical-align:bottom">.
+    * init_graph_state = []. A list of integer tuples that specify the qubit pairs for performing a controlled-Z gate on, to generate an initial graph state, after first initializing all along the +x axis. If 'init_graph_state' is used, all other initialization parameters should be left empty. The qubit pairs will represent j and k in the graph state formula $\left|\psi_0\rangle= \prod_{(j,k)\in V}{CZ}[j,k] \prod_i \right|+x_i\rangle.$
 	* init_pauli_state = "": This initialization parameter is deprecated, use 'init_product_state' instead.
     * load_files_prefix = "" (str): The prefix of files as previously saved using the simulator, which the initial state has to be loaded from. An empty string indicates that the initial state is not loaded. If 'load_files_prefix' is used, all other initialization parameters should be left empty. See the parameter 'b_save_final_state' for more details on the saved files.
 * Additional operations:
     * apply_gates = []. A list of tuples of the form `(time: float, gate: str, q0: int, q1: Optional[int])` that specify one-qubit or two-qubit instantaneous gates that will be applied at the specified times, on the indicated qubits. At each time, the order of application of the gates is according to their order in the list. The supported gate strings are the three Paulis ("X", "Y", "Z"), Hadamard ("H"), Sqrt-X ("SX") Controlled-X (CNOT, "CX") and Controlled-Z ("CZ").
     * b_apply_gate_compression = True (bool): Whether to do a state truncation with parameter `cut_off_rho` after each two-qubit gate.
 * Lattice specification:
-    * l_x = 0 (float): The length of the lattice along the x dimension. In case of value 0, the number of qubits <img src="https://render.githubusercontent.com/render/math?math=N" style="vertical-align:bottom"> is used, and parameter l_y must be 1.
+    * l_x = 0 (float): The length of the lattice along the x dimension. In case of value 0, the number of qubits N is used, and parameter l_y must be 1.
     * l_y = 1 (float): The length of the lattice along the y dimension.
     * b_periodic_x = False (bool): Whether periodic boundary conditions are applied along the x dimension. If True, then l_y must be 1. If False, open boundary conditions are used along the x dimension.
     * b_periodic_y = False (bool): Whether periodic boundary conditions are applied along the y dimension. If False, open boundary conditions are used along the y dimension.
 * Numerical simulation control:
     * trotter_order = 4 (int): Trotter approximation order, Possible values are 2, 3, 4.
     * max_dim_rho = 400 (int): Maximum bond dimension for density matrices.
     * cut_off_rho = 1e-16 (float): Maximum truncation error (discarded Schmidt weight) for density matrices. The actual truncation is done using the most severe condition between cut_off_rho and max_dim_rho.
-    * b_force_rho_trace = True (bool): Whether to force the density matrix trace to one by substituting <img src="https://render.githubusercontent.com/render/math?math=\rho \to\rho/ {\rm tr}\{\rho\}" style="vertical-align:bottom"> at every time step, compensating for some finite-step errors.
-    * force_rho_hermitian_step = 4 (int): Determines every how many evolution time steps (<img src="https://render.githubusercontent.com/render/math?math=\tau" style="vertical-align:bottom">), to substitute <img src="https://render.githubusercontent.com/render/math?math=\rho \to (\rho %2b \rho^\dagger)/2" style="vertical-align:bottom">. This may reduce some errors, but is computationally expensive.
+    * b_force_rho_trace = True (bool): Whether to force the density matrix trace to one by substituting $\rho \to\rho/ {\rm tr}\{\rho\}$ at every time step, compensating for some finite-step errors.
+    * force_rho_hermitian_step = 4 (int): Determines every how many evolution time steps ($\tau$), to substitute $\rho \to (\rho + \rho^\dagger)/2$. This may reduce some errors, but is computationally expensive.
     * b_initial_rho_compression = True (bool): Whether a density matrix that is loaded from a previously saved state, should be re-gauged and compressed. Has no effect if the initial state is not loaded from a previously saved state.
 * Observables and output:
     * 1q_indices = [] (list[int]): A list of integers that specify the qubits for which single-qubit observables will be calculated. In the case of an empty list, single-qubit observables will be calculated for all qubits.
     * 1q_components = ['Z'] (list[str]): A list of strings that specify the Pauli observables to compute for all qubits given in parameter "1q_indices". The allowed strings in the list are "x", "y", or "z" (lower or upper case). The observables results are saved using a file name ending with "obs-1q.dat".
     * 2q_indices = [] (list[tuple(int)]): A list of integer tuples that specify the qubit pairs for calculating two-qubit expectation values. In the case of an empty list, two-qubit expectation values will be calculated for all qubit pairs.
     * 2q_components = ['ZZ'] (list[str]): A list of strings that specify the two-qubit Pauli observables to compute for all qubit pairs given in parameter "2q_indices". The allowed strings in the list are one of "xx", "yy", "zz", "xy", "xz", "yz", "yx", "zx", "zy" (lower or upper case). The observables results are saved using a file name ending with ".obs-2q.dat".
-    * output_step = 1 (int): How often (in integer steps of time <img src="https://render.githubusercontent.com/render/math?math=\tau" style="vertical-align:bottom">) the observables are computed. In case of the value 0, no observables will be computed.
+    * output_step = 1 (int): How often (in integer steps of time $\tau$) the observables are computed. In case of the value 0, no observables will be computed.
     * b_save_final_state = False (bool): Whether to save the final state to files (Three binary files will be saved, whose names will have the prefix defined in "output_files_prefix").
     * b_quiet = False (bool): Whether to avoid writing the console output at every time step (while the output at every time step is still written to the log file). The initialization step output and the final output information are always written to the console (together with the log file).