We use the same parameter set as used in Ref. [1]}: \(B_z = 8 \eta \) , \(\Lambda = 100 \eta \) , \(\omega _0 = 2\eta \) , and \(\beta = 0.5/ \eta \)  for various \(\epsilon ^2\)  within the range \( 0.02 \eta \le \epsilon ^2 \le 2 \eta \) . The system is initially in the state with all spins aligned against \(B_z\) . As we carried out in Ref. [1]}, we ignore the Lamb shift term. With these conditions, we numerically obtain the steady-state density matrix \(\rho ^{\rm ss}\)  by directly solving the ULE for the spin chain with length \(N = 6\) . As expected, the obtained steady state deviates from the simple Gibbs state of the system Hamiltonian as shown in Fig. REF (a). Since the Lamb shift term is ignored, \(\rho _{nn}^{(2)}\)  term should vanish according to the analytic continuation method. Therefore, the next order of the diagonal term of the steady-state reduced density matrix for this ULE is \(O(\epsilon ^4)\) , i.e., \(\rho _{ nn}^{\rm ss} = (\rho _{\rm G})_{nn} + O(\epsilon ^4)\) . To verify this numerically, we make a plot of \(\Delta \rho _{nn} \equiv [\rho _{nn}^{\rm ss} - (\rho _{\rm G})_{nn}]/ (\rho _{\rm G})_{nn}\)  as a function of \(\epsilon ^2\)  for various \(n\)  and check whether \(\Delta \rho _{nn} \sim \epsilon ^4\)  as expected from the analytic continuation. The result is presented in Fig. REF (b). As seen in the figure, \(\Delta \rho _{nn} \)  shows \(\epsilon ^4\)  behavior for small \(\epsilon \) . This clear demonstrates that the analytic calculation provides the correct answer.
<FIGURE>