Coarse Graining Details

By "coarse graining", we mean fitting the spectral function using a finite number of broad bosonic modes. It means that we limit the number of bath modes to some finite number, but introduce broadenings \(\omega_i\rightarrow\omega_i+\textrm{i}\gamma/2\). A similar approach has been taken in many classical numerical methods, such as in the Hierachy of Equations of Motion (HEOM), or various methods of introducing auxiliary modes with Lindblad decay.

Broad Boson Modes

Our goal is to construct a bath that reproduces the desired spectral function as well as possible, but includes only a finite number of boson modes. The approach is to broaden the bosonic modes so that they effectively represent a continuous bath. Such a broadening can be established theoretically by introducing a damping Lindbladian:

\[ \dot \rho = \textrm{i}[\rho, H] + \gamma\left( a\rho a^\dagger -\frac{1}{2}a^\dagger a\rho -\frac{1}{2}\rho a^\dagger a \right) \nobreakspace . \]

On the Hamiltonian level, the equivalent is to introduce a bilinear coupling between a boson mode and a "flat" environment (bath of the bath). This coupling leads to a broadening of the spectral peaks,

\[ \begin{align} \langle a(t)a^\dagger\rangle &= e^{-\textrm{i}\omega_0 t - \gamma \vert t \vert/2} \\ \int_{-\infty}^{\infty} e^{\textrm{i}\omega t}\langle a(t)a^\dagger\rangle dt &= 2\pi \frac{1}{\pi}\frac{\gamma/2}{(\gamma/2)^2+(\omega-\omega_0)^2} \end{align} \nobreakspace . \]

Coupling such modes to a system operator such as \(\sigma_z/2\) via a coupling coefficient \(g\) gives the corresponding spectral function peak:

\[ S(\omega) = 2\pi g^2 \frac{1}{\pi}\frac{\gamma/2}{(\gamma/2)^2+(\omega-\omega_0)^2} \nobreakspace . \]

Broad Boson Modes Represented by Noisy Spins

The coarse-grained bath is bosonic. We aim to represent this using a spin bath. Various efficient "digital" (e.g. binary or unary) mappings between qubit states and boson states have been proposed in recent literature. However, such encodings do not map the decay of (arbitrary) bath qubits to single-boson annihilation,which is the key correspondence in our approach.

Instead, the spin-boson simulator directly replaces the coarse-grained bath operators with spin-bath operators, in the simplest case by doing the following:

\[ \begin{align} g(a + a^\dagger) &\rightarrow g \sigma_x \\ \omega a^\dagger a &\rightarrow \frac{\omega}{2}\sigma_z \\ \sqrt{\gamma} a &\rightarrow \sqrt{\gamma}\sigma^- \end{align} \]

This can work, since spins can (collectively) behave like a bosonic field. The condition for this to be true boils down to the requirement that the coupling-operator statistics to be Gaussian, up to some relevant order. This is clearly valid for a spin-bath with a fast memory decay compared to the system-bath coupling, \(\gamma\gg g\), since here only two-operator correlation functions matter. A less trivial example is a group of identical spins that can collectively behave like one bosonic mode with long memory time. Such a system is explicitly studied in this example.

More generally, we do not have to assume that the spins are identical. A central condition for the bath to be bosonic is that individual bath spins are only weakly excited,

\[ \left< \sigma^+_i \sigma^-_i \right> _{\textrm{full solution}} \ll 1 \quad i\in \textrm{bath} \nobreakspace . \]

Note that the collective excitation number, \( \sum_{i\in \textrm{bath}}\left< \sigma^+_i \sigma^-_i \right> _\textrm{full solution}\), can still be fairly large. This relation should be valid for all times during the simulation. It can be probed by bath-qubit measurements.

If this demand is not respected, the "bosonicity" of the bath can be improved by representing a problematic mode by \(N \) spins, all of them with a down-scaled coupling:

\[ g \rightarrow \frac{g}{\sqrt{N}} \nobreakspace . \]

Downscaling guarantees that the the total spectral function remains unchanged. We then use:

\[ \begin{align} g(a + a^\dagger) &\rightarrow \sum_{i=1}^N \frac{g}{\sqrt{N}} \sigma^x_i \\ \omega a^\dagger a &\rightarrow \sum_{i=1}^N \frac{\omega}{2}\sigma^z_i \\ \sqrt{\gamma}a & \rightarrow \sum_{i=1}^N \sqrt{\gamma}\sigma^-_i \end{align} \] It should be noted that it can also be enough to just increase the total number of broad modes in the spectral-function fitting range (keeping \(N=1\)).

In coarse-graining theory, the broadening of the bosonic modes is always due to damping. However, the corresponding broadening of the spin modes can also be partly due to dephasing. The correspondence between the broadenings is straightforward:

\[ \gamma = \gamma_\textrm{damping} + \gamma_\textrm{dephasing} \nobreakspace . \]

Noiseless System

In the fitting, our target function is the frequency-discretized multidimensional spectral function:

\[ S_{ij}(\omega_m) \equiv S[i, j, m] , . \] It is spin-exchange symmetric, so \(S_{ij}(\omega_m)= S_{ji}(\omega_m)\). More generally, it is Hermitian, as we are dealing with complex couplings.

Our fitting finds an approximation of this function by replacing the original bath with \(n_\textrm{lorentzians}\) broad modes. Each broad mode has a central frequency \(\omega_k\), coupling to the \(i^\text{th}\) spin \(g_{ik}\), and width \(\gamma_k\). The (multidimensional) spectral function that follows from this is:

\[ S_{ij}^\textrm{cg}(\omega) = 2\pi\sum_{k=1}^{n_\textrm{lorentzians}} g_{ik}g_{jk}^*\frac{1}{\pi}\frac{\gamma_k/2}{(\gamma_k/2)^2+(\omega-\omega_k)^2} \nobreakspace . \] The current fitting routine optimizes the following cost function:

\[ C= \sum_{i\leq j} \left \lbrace\sum_k \left [S_{ij}^\textrm{cg}(\omega_k) - S_{ij}(\omega_k) \right]^2 \right \rbrace \nobreakspace , \]

by searching for the optimal values of \(\omega_k\), \(g_{ik}\), and \(\gamma_k\), starting from some initial guess. The general fitting problem then has \(n_\textrm{lorentzians}\) frequencies, \(n_\textrm{spins} \times n_\textrm{lorentzians}\) couplings, and \(n_\textrm{lorentzians}\) widths.

The fitting can also be done using the same broadening for all modes, \(\gamma_m=\gamma\). This may be a more realistic choice when mapping to hardware properties.

Note that if \(n_\textrm{spins}=1\), we are dealing with fitting only one function, \(S_{ij}(\omega_m) \equiv S(\omega_m)\). Furthermore, if \(n_\textrm{spins}>1\), but cross-correlations are very small and diagonal entries (\(i=j\)) of the spectrum are almost identical, it is recommended to perform a comparison with the result of fitting only the corresponding single-spin systems.

Open System

We also allow the user to add system-qubit noise as a background to the spectral function and perform the fitting taking this into account. This is motivated by the possibility to add system damping exactly to fermion-transport models. We then introduce a system rate that is an average of bath rates, times some given factor \(f\),

\[ \gamma_\textrm{system} \equiv f \langle \gamma_i\rangle , , \]

and change the spectral function as:

\[ S_{ij}^\textrm{cg}(\omega) = 2\pi\sum_{k=1}^{n_\textrm{lorentzians}} g_{ik}g_{jk}^*\frac{1}{\pi}\frac{\gamma_k}{(\gamma_k/2)^2+(\omega-\omega_k)^2} + \delta_{ij} \gamma_\textrm{system} \nobreakspace . \]

The noise of different system qubits is assumed to be uncorrelated, hence the multiplication by \(\delta_{ij}\). The correct factor \(f\) to be used depends on the model considered (as well as on hardware properties). Note that the fitting problem remains essentially the same.

Spectral Function to Coupling

Here, we assume that a spectral function is given as an input and we want to write down a Hamiltonian that reproduces it: we need to solve the system-bath couplings from the spectral function. The spectral function has been discretized to frequencies \(\omega_m\), and the corresponding data values are \(S_{ij}(\omega_m)\equiv S[i, j, m]\). In the corresponding Hamiltonian description, the frequencies of the bath modes will be \(\omega_m\). It should be noted that if the given spectral function is unphysical, the software will find the closest physical solution (according to the Frobenius norm).

In the case of a single-spin system, the couplings are defined by the areas:

\[ S(\omega_m) d\omega_m = |g_m|^2 \nobreakspace , \]

where \(d\omega_m\) is the width of the frequency interval corresponding to mode \(\omega_m\). The couplings are then solved to be \(g_m=\pm\sqrt{S(\omega_m) d\omega_m}\). The chosen sign is not important here.

In the case of \(n\) system spins, the couplings are constructed similarly with the help of linear algebra. We introduce \(n\) bath-submodes per frequency and solve:

\[ S[i, j, m] d\omega_m = \sum_{s=1}^{n} g[i, s, m] g*[j, s, m] \nobreakspace . \]

Here \(g[i, s, m]\) is the coupling from spin \(i\) to bath submode \(s\) at frequency \(\omega_m\). The spectral function is Hermitian with respect to the changing \(i\) and \(j\). There are multiple solutions to this equation, some of which have unreasonable and complicated cross-couplings. To find a reasonable solution, we can require:

\[ g[i, s, m] = 0 \qquad \textrm{for } i < s \ . \] In this form, the obtained cross-couplings are zero if there are no spectral-function cross-correlations. The solution can then be found by Cholesky decomposition:

\[ S[:, :, m] = L L^\dagger \ , \] where \(L\) is the solution for the coupling matrix at frequency \(\omega_m\).