  We present an essentially exact numerical method for modelling flux qubit
chains subject to charge and flux noise. We define an essentially exact method
as one that introduces errors that are completely controlled such that they can
be made arbitrarily small by tuning the simulation parameters. The method
adopts the quasi-adiabatic path integral formalism to express the system's
reduced density matrix as a time-discretized path integral, comprising a series
of influence functionals that encode the non-Markovian dynamics of the system.
We present a detailed derivation of the path integral expression for the
system's reduced density matrix and describe in detail the tensor network
algorithm used to evaluate the path integral expression. We have implemented
our method in an open-sourced Python library called "spinbosonchain". When
appropriate, we draw connections between concepts covered in this manuscript
and the library's code.
