diff --git a/.github/scripts/fix-timestamps b/.github/scripts/fix-timestamps
new file mode 100644
index 0000000..d6dd630
--- /dev/null
+++ b/.github/scripts/fix-timestamps
@@ -0,0 +1,9 @@
+#!/usr/bin/env sh
+
+# This script fixes the timestamps of all files in the current git repository
+
+git ls-tree -r --name-only HEAD | while read filename; do
+  unixtime=$(git log -1 --format="%at" -- "${filename}")
+  touchtime=$(date -d @$unixtime +'%Y%m%d%H%M.%S')
+  touch -t ${touchtime} "${filename}"
+done
\ No newline at end of file
diff --git a/.github/workflows/preview.yml b/.github/workflows/preview.yml
index 5dbd3e7..19db5ff 100644
--- a/.github/workflows/preview.yml
+++ b/.github/workflows/preview.yml
@@ -29,14 +29,10 @@ jobs:
       - uses: julia-actions/julia-buildpkg@v1
         with:
           precompile: true
-      - uses: cardinalby/export-env-action@v2
-        with:
-          envFile: '_environment'    
-          expand: 'true'
+      - name: Fix Timestamps
+        run: bash .github/scripts/fix-timestamps
       - name: Quarto Render
         uses: quarto-dev/quarto-actions/render@v2
-        with:
-          to: html
       - name: Deploy PR Preview
         uses: rossjrw/pr-preview-action@v1.6.1
         with:
diff --git a/.github/workflows/publish.yml b/.github/workflows/publish.yml
index 97989cf..a3208a8 100644
--- a/.github/workflows/publish.yml
+++ b/.github/workflows/publish.yml
@@ -7,7 +7,7 @@ on: # either one of the following three cases
     - cron: '0 0 * * 0' # weekly (every Sunday)
 
   push:
-    branches: main
+    branches: [main]
 
 jobs:
   build:
@@ -26,10 +26,8 @@ jobs:
       - uses: julia-actions/julia-buildpkg@v1
         with:
           precompile: true
-      - uses: cardinalby/export-env-action@v2
-        with:
-          envFile: '_environment'    
-          expand: 'true'
+      - name: Fix Timestamps
+        run: bash .github/scripts/fix-timestamps
       - name: Quarto Render
         uses: quarto-dev/quarto-actions/render@v2
         with:
diff --git a/HierarchicalEOM.jl/SIAM.qmd b/HierarchicalEOM.jl/SIAM.qmd
index 70d0b14..df37e85 100644
--- a/HierarchicalEOM.jl/SIAM.qmd
+++ b/HierarchicalEOM.jl/SIAM.qmd
@@ -1,7 +1,7 @@
 ---
 title: "Single-impurity Anderson model"
 author: Yi-Te Huang
-date: 2025-05-22  # last update (keep this comment as a reminder)
+date: last-modified
 ---
 
 ## Introduction
@@ -27,6 +27,8 @@ Now, we need to build the system Hamiltonian and initial state with the package
 ```{julia}
 using HierarchicalEOM  # this automatically loads `QuantumToolbox`
 using CairoMakie       # for plotting results
+
+CairoMakie.activate!(type = "svg")
 ```
 
 ```{julia}
@@ -88,7 +90,7 @@ dos = DensityOfStates(M_odd, ados_s, d_up, ωlist)
 plot the results
 
 ```{julia}
-fig = Figure(size = (500, 350))
+fig = Figure()
 ax = Axis(fig[1, 1], xlabel = L"\omega")
 lines!(ax, ωlist, dos)
 
diff --git a/HierarchicalEOM.jl/cavityQED.qmd b/HierarchicalEOM.jl/cavityQED.qmd
index e4a191b..478c75a 100644
--- a/HierarchicalEOM.jl/cavityQED.qmd
+++ b/HierarchicalEOM.jl/cavityQED.qmd
@@ -1,7 +1,7 @@
 ---
 title: "Cavity QED system"
 author: Shen-Liang Yang, Yi-Te Huang
-date: 2025-05-22  # last update (keep this comment as a reminder)
+date: last-modified
 ---
 
 ## Introduction
@@ -41,6 +41,8 @@ Now, we need to build the system Hamiltonian and initial state with the package
 ```{julia}
 using HierarchicalEOM  # this automatically loads `QuantumToolbox`
 using CairoMakie       # for plotting results
+
+CairoMakie.activate!(type = "svg")
 ```
 
 ```{julia}
@@ -96,7 +98,7 @@ Ct = correlation_function(Bath, tlist_test);
 Ct2 = correlation_function(Bath_test, tlist_test)
 
 # plot
-fig = Figure(size = (500, 350))
+fig = Figure()
 ax = Axis(fig[1, 1], xlabel = L"t", ylabel = L"C(t)")
 lines!(ax, tlist_test, real(Ct2), label = L"$N=1000$ (real part)", linestyle = :solid)
 lines!(ax, tlist_test, real(Ct),  label = L"$N=20$ (real part)", linestyle = :dash)
@@ -154,7 +156,7 @@ np_steady_H = expect(a' * a, steady_H)
 plot results
 
 ```{julia}
-fig = Figure(size = (600, 350))
+fig = Figure()
 
 ax1 = Axis(fig[1, 1], xlabel = L"t")
 lines!(ax1, t_list, σz_evo_H, label = L"\langle \sigma_z \rangle", linestyle = :solid)
@@ -176,7 +178,7 @@ fig
 psd_H = PowerSpectrum(M_Heom, steady_H, a, ω_list)
 
 # plot
-fig = Figure(size = (500, 350))
+fig = Figure()
 ax = Axis(fig[1, 1], xlabel = L"\omega")
 lines!(ax, ω_list, psd_H)
 
@@ -223,7 +225,7 @@ np_steady_M = expect(a' * a, steady_M);
 plot results
 
 ```{julia}
-fig = Figure(size = (600, 350))
+fig = Figure()
 
 ax1 = Axis(fig[1, 1], xlabel = L"t")
 lines!(ax1, t_list, σz_evo_M, label = L"\langle \sigma_z \rangle", linestyle = :solid)
diff --git a/HierarchicalEOM.jl/dynamical_decoupling.qmd b/HierarchicalEOM.jl/dynamical_decoupling.qmd
index 30e89b5..6180d34 100644
--- a/HierarchicalEOM.jl/dynamical_decoupling.qmd
+++ b/HierarchicalEOM.jl/dynamical_decoupling.qmd
@@ -1,7 +1,7 @@
 ---
 title: "Driven systems and dynamical decoupling"
 author: Yi-Te Huang
-date: 2025-05-22  # last update (keep this comment as a reminder)
+date: last-modified
 ---
 
 Inspirations taken from an example in QuTiP-BoFiN article [@QuTiP-BoFiN2023].
@@ -37,6 +37,8 @@ Now, we need to build the system Hamiltonian and initial state with the package
 ```{julia}
 using HierarchicalEOM  # this automatically loads `QuantumToolbox`
 using CairoMakie       # for plotting results
+
+CairoMakie.activate!(type = "svg")
 ```
 
 ```{julia}
@@ -91,7 +93,7 @@ fastTuple = (V = amp_fast, Δ = delay)
 slowTuple = (V = amp_slow, Δ = delay)
 
 # plot
-fig = Figure(size = (600, 350))
+fig = Figure()
 ax = Axis(fig[1, 1], xlabel = L"t")
 lines!(ax, tlist, [pulse(fastTuple, t) for t in tlist], label = "Fast Pulse", linestyle = :solid)
 lines!(ax, tlist, [pulse(slowTuple, t) for t in tlist], label = "Slow Pulse", linestyle = :dash)
@@ -152,7 +154,7 @@ slowPulseSol = HEOMsolve(M, ψ0, tlist; e_ops = [ρ01], H_t = H_D, params = slow
 ## Plot the coherence
 
 ```{julia}
-fig = Figure(size = (600, 350))
+fig = Figure()
 ax = Axis(fig[1, 1], xlabel = L"t", ylabel = L"\rho_{01}")
 lines!(ax, tlist, real(fastPulseSol.expect[1, :]), label = "Fast Pulse", linestyle = :solid)
 lines!(ax, tlist, real(slowPulseSol.expect[1, :]), label = "Slow Pulse", linestyle = :dot)
diff --git a/HierarchicalEOM.jl/electronic_current.qmd b/HierarchicalEOM.jl/electronic_current.qmd
index 08ba820..2bf87ab 100644
--- a/HierarchicalEOM.jl/electronic_current.qmd
+++ b/HierarchicalEOM.jl/electronic_current.qmd
@@ -1,7 +1,7 @@
 ---
 title: "Electronic Current"
 author: Yi-Te Huang
-date: 2025-05-22  # last update (keep this comment as a reminder)
+date: last-modified
 ---
 
 Inspirations taken from [qutip documentation](https://qutip.org/documentation.html).
@@ -28,6 +28,8 @@ Now, we need to build the system Hamiltonian and initial state with the package
 ```{julia}
 using HierarchicalEOM  # this automatically loads `QuantumToolbox`
 using CairoMakie       # for plotting results
+
+CairoMakie.activate!(type = "svg")
 ```
 
 ```{julia}
@@ -149,7 +151,7 @@ end
 plot the result
 
 ```{julia}
-fig = Figure(size = (500, 350))
+fig = Figure()
 ax = Axis(fig[1, 1], xlabel = "time", ylabel = "Current")
 lines!(ax, tlist, Ie_L, label = "Bath L", color = :blue, linestyle = :solid)
 lines!(ax, tlist, Ie_R, label = "Bath R", color = :red, linestyle = :solid)
diff --git a/QuantumToolbox.jl/time_evolution/Dicke.qmd b/QuantumToolbox.jl/time_evolution/Dicke.qmd
index b1c29fb..3f7785c 100644
--- a/QuantumToolbox.jl/time_evolution/Dicke.qmd
+++ b/QuantumToolbox.jl/time_evolution/Dicke.qmd
@@ -1,7 +1,7 @@
 ---
 title: "The Dicke Model"
 author: Li-Xun Cai
-date: 2025-05-22  # last update (keep this comment as a reminder)
+date: last-modified
 ---
 
 Inspirations taken from [this QuTiP tutorial](https://nbviewer.org/urls/qutip.org/qutip-tutorials/tutorials-v5/lectures/Lecture-3A-Dicke-model.ipynb) by J. R. Johansson.
@@ -48,6 +48,8 @@ This formulation reduces complexity, as it allows us to work on a collective bas
 ```{julia}
 using QuantumToolbox
 using CairoMakie
+
+CairoMakie.activate!(type = "svg")
 ```
 
 ```{julia}
@@ -105,7 +107,7 @@ Jzvec = expect(Jz0, ψGs);
 ```
 
 ```{julia}
-fig = Figure(size = (800, 300))
+fig = Figure(size=(900, 350))
 axn = Axis(
  fig[1,1],
  xlabel = "interaction strength",
@@ -119,7 +121,7 @@ axJz = Axis(
 ylims!(-N0/2, N0/2)
 lines!(axn, gs, real(nvec))
 lines!(axJz, gs, real(Jzvec))
-display(fig);
+fig
 ```
 The expectation value of photon number and $\hat{J}_z$ showed a sudden increment around $g_c$. 
 
@@ -140,11 +142,12 @@ for (hpos, idx) in enumerate(cases)
     
     # plot fock distribution
     _, ax2 = plot_fock_distribution(ρcav, location = fig[2,hpos])
+
+    ax2.xticks = (0:2:size(ρcav,1)-1, string.(0:2:size(ρcav,1)-1))
     
     if hpos != 1
-        ax.xlabelvisible, ax.ylabelvisible, ax2.ylabelvisible = fill(false, 3)
-        ax.xticksvisible, ax.yticksvisible, ax2.yticksvisible = fill(false, 3)
-        ax2.yticks = (0:0.5:1, ["" for i in 0:0.5:1])
+        hideydecorations!(ax, ticks=false)
+        hideydecorations!(ax2, ticks=false)
         if hpos == 5 # Add colorbar with the last returned heatmap (_hm) 
             Colorbar(fig[1,6], hm)
         end
@@ -158,7 +161,7 @@ lines!(ax3, gs, real(nvec), color=:teal)
 ax3.xlabelsize, ax3.ylabelsize = 20, 20
 vlines!(ax3, gs[cases], color=:orange, linestyle = :dash, linewidth = 4)
 
-display(fig);
+fig
 ```
 As $g$ increases, the cavity ground state's wigner function plot looks more coherent than a thermal state.
 
@@ -174,7 +177,7 @@ end;
 ```
 
 ```{julia}
-fig = Figure(size=(800, 400))
+fig = Figure(size=(900, 400))
 ax = Axis(fig[1,1])
 ax.xlabel = "coupling strength"
 ax.ylabel = "mutual entropy"
@@ -184,7 +187,7 @@ for (idx, slist) in enumerate(slists)
 end
 
 Legend(fig[1,2], ax, label = "number of atoms")
-display(fig);
+fig
 ```
 We further consult mutual entropy between the cavity and the spins as a measure of their correlation; the result showed that as the number of atoms $N$ increases, the peak of mutual entropy moves closer to $g_c$.
 
diff --git a/QuantumToolbox.jl/time_evolution/adiabatic.qmd b/QuantumToolbox.jl/time_evolution/adiabatic.qmd
index 86bbd64..9fb747b 100644
--- a/QuantumToolbox.jl/time_evolution/adiabatic.qmd
+++ b/QuantumToolbox.jl/time_evolution/adiabatic.qmd
@@ -1,7 +1,7 @@
 ---
 title: "Adiabatic sweep (with `QuantumObjectEvolution`)"
 author: Li-Xun Cai
-date: 2025-05-22  # last update (keep this comment as a reminder)
+date: last-modified
 ---
 
 Inspirations taken from [the QuTiP tutorial](https://nbviewer.org/urls/qutip.org/qutip-tutorials/tutorials-v5/lectures/Lecture-8-Adiabatic-quantum-computing.ipynb) by J. R. Johansson.
@@ -40,6 +40,8 @@ where the parameter $T$ determines how slow the Hamiltonian changes in time.
 ```{julia}
 using QuantumToolbox
 using CairoMakie
+
+CairoMakie.activate!(type = "svg")
 ```
 
 ```{julia}
@@ -104,7 +106,7 @@ end
 ```
 
 ```{julia}
-fig = Figure(size=(800, 400))
+fig = Figure(size=(900, 400))
 ax = Axis(fig[1,1], xticks = (0:0.25:1, ["$(t)T" for t in 0:0.25:1]))
 
 for idx in 1:20 # only check for the lowest 20 eigenvalues
@@ -112,7 +114,7 @@ for idx in 1:20 # only check for the lowest 20 eigenvalues
     lines!(ax, range(0,1,length(tlist)), eigs[:,idx], label = string(idx), color = color)
 end
 
-display(fig)
+fig
 ```
 
 The plot shows that the gap is nonvanishing and thus validates the evolution. So we proceed to check the expectation value dynamics of the final Hamiltonian and the fidelity dynamics to the truthful ground state throughout the evolution.
@@ -127,7 +129,7 @@ end;
 ```
 
 ```{julia}
-fig = Figure(size=(800, 400))
+fig = Figure(size=(900, 400))
 axs = Axis.([fig[1,1], fig[1,2]])
 axs[1].title = L"\langle H_f \rangle"
 axs[1].xticks = (0:0.25:1, ["$(t)T" for t in 0:0.25:1])
@@ -143,7 +145,7 @@ end
 
 Legend(fig[1,3], axs[1], L"T")
 
-display(fig)
+fig
 ```
 As the plot showed, the fidelity between the prepared final state and the truthful ground state reaches 1 as the total evolution time $T$ increases, showcasing the requirement of the adiabatic theorem that the change has to be gradual. 
 
diff --git a/QuantumToolbox.jl/time_evolution/fluorescence.qmd b/QuantumToolbox.jl/time_evolution/fluorescence.qmd
index 2c30211..5a5bb70 100644
--- a/QuantumToolbox.jl/time_evolution/fluorescence.qmd
+++ b/QuantumToolbox.jl/time_evolution/fluorescence.qmd
@@ -1,7 +1,7 @@
 ---
 title: "Resonance fluorescence"
 author: Li-Xun Cai
-date: 2025-05-22 # last update (keep this comment as a reminder)
+date: last-modified
 ---
 
 Inspirations taken from [this QuTiP tutorial](https://nbviewer.org/urls/qutip.org/qutip-tutorials/tutorials-v5/lectures/Lecture-13-Resonance-flourescence.ipynb) by J. R. Johansson.
@@ -81,7 +81,9 @@ where
 
 ```{julia}
 using QuantumToolbox
-import CairoMakie: Figure, Axis, @L_str, lines!, axislegend, display, ylims!, xlims!
+using CairoMakie
+
+CairoMakie.activate!(type = "svg")
 ```
 
 ```{julia}
@@ -116,7 +118,7 @@ e_ops = [
 ]
 ψ0 = e_ket # set initial state being purely excited to better observe the radiative behaviour
 L = liouvillian_spec(Ω, γ0, KT)
-print(L)
+L
 ```  
 
 We already generate the Liouvillian with `c_ops` included above. We don't need to specify the `c_ops` again in [`mesolve`](https://qutip.org/QuantumToolbox.jl/stable/resources/api#QuantumToolbox.mesolve):  
@@ -128,7 +130,7 @@ sol = mesolve(L, ψ0, tlist, nothing, e_ops = e_ops)
 By observing the expectation values of the Pauli operators, we see that the Bloch vector $(\langle \hat{\sigma}_x \rangle, \langle \hat{\sigma}_y \rangle, \langle \hat{\sigma}_z \rangle)$ becomes shorter over time, which is consistent with the dissipative behaviour. Also, the population of the excited state $|e\rangle$ has an oscillation amplitude decaying over time.  
 ```{julia}
 expect = real.(sol.expect)
-fig1 = Figure(size = (600,300))
+fig1 = Figure(size = (900,400))
 ax11 = Axis(
     fig1[1,1]
 )
@@ -146,7 +148,7 @@ lines!(ax12, tlist, expect[6,:], label = L"P_e")
 axislegend(ax12)
 ylims!(ax12, (0,1))
 
-display(fig1);
+fig1
 ```
 
 Further, we check the effect of different values of the damping rate. Note that despite these dissipation rates looked enormous at first glance, it is still to the order of the field strength and therefore considered dissipative for the system in the lab frame.
@@ -162,7 +164,7 @@ end
 
 The expectation values dynamics of $\hat{\sigma}^{+}$ and $\hat{\sigma}^{-}$ shows the driving-field-induced dipole moment of the atom oscillates and persists.
 ```{julia}
-fig2 = Figure(size = (600,300))
+fig2 = Figure(size = (900,400))
 ax2 = Axis(
     fig2[1,1],
     xlabel = L"time $[1/\Omega]$",
@@ -174,11 +176,11 @@ for (γ0, expect) in results
 end
 
 axislegend(ax2)
-display(fig2);
+fig2
 ```
 
 ```{julia}
-fig3 = Figure(size = (600,300))
+fig3 = Figure(size = (900,400))
 ax3 = Axis(
     fig3[1,1],
     xlabel = L"time $[1/\Omega]$",
@@ -189,12 +191,12 @@ for (γ0, expect) in results
     lines!(ax3, tlist, imag(expect[4,:]), label = "γ0 = $γ0")
 end
 axislegend(ax3)
-display(fig3);
+fig3
 ```
 
 We now move to the analysis of the correlation function $C(\tau) = \langle \hat{\sigma}^{+}(\tau) \hat{\sigma}^{-}(0)\rangle$, which describes the radiative behaviour of the atom towards its surrounding environment. Using [`correlation_2op_1t`](https://qutip.org/QuantumToolbox.jl/stable/resources/api#QuantumToolbox.correlation_2op_1t), we can obtain the correlation function as a function of $\tau$ and use [`spectrum_correlation_fft`](https://qutip.org/QuantumToolbox.jl/stable/resources/api#QuantumToolbox.spectrum_correlation_fft) to obtain the corresponding Fourier transform.
 ```{julia}
-fig4 = Figure(size = (600,300))
+fig4 = Figure(size = (900,400))
 ax41 = Axis(
     fig4[1,1],
     xlabel = L"\tau",
@@ -219,7 +221,7 @@ end
 xlims!(ax42, (-2,2))
 axislegend(ax41)
 axislegend(ax42)
-display(fig4);
+fig4
 ```
 In the above plots, one finds that the correlation functions decay faster with higher dissipation rate, and therefore the lower spectral peaks. On the other hand, the higher spectral peaks means the radiation is brighter in terms of intensity.
 
diff --git a/QuantumToolbox.jl/time_evolution/kerr.qmd b/QuantumToolbox.jl/time_evolution/kerr.qmd
index 30d4a96..54066cd 100644
--- a/QuantumToolbox.jl/time_evolution/kerr.qmd
+++ b/QuantumToolbox.jl/time_evolution/kerr.qmd
@@ -1,7 +1,7 @@
 ---
 title: "Kerr nonlinearities"
 author: Li-Xun Cai
-date: 2025-05-22  # last update (keep this comment as a reminder)
+date: last-modified
 ---
 
 Inspirations taken from [this QuTiP tutorial](https://nbviewer.org/urls/qutip.org/qutip-tutorials/tutorials-v5/lectures/Lecture-14-Kerr-nonlinearities.ipynb) by J. R. Johansson.
@@ -54,6 +54,8 @@ where $\chi$ again absorbed the coefficients.
 ```{julia}
 using QuantumToolbox
 using CairoMakie
+
+CairoMakie.activate!(type = "svg")
 ```
 
 We begin by defining functions for visualization:
@@ -73,7 +75,7 @@ function plot_variance(op, tlist, states)
 end
 
 function plot_Fock_dist(tlist, states)
-    fig = Figure()
+    fig = Figure(size = (900, 300))
     ax = Axis(
         fig[1,1],
         xlabel = L"N",
@@ -168,7 +170,7 @@ idx = searchsortedfirst(tlist, π/χ) - 1
 fig3, ax3, hm3 = plot_wigner(result.states[idx])
 ax3.title = "χt = $(tlist[idx])"
 Colorbar(fig3[1,2], hm3)
-display(fig3)
+fig3
 ```
 As the plot revealed, the state at $\chi t = \pi$ is in fact the [**cat state**](https://en.wikipedia.org/wiki/Cat_state#Cat_states_in_single_modes) for a single mode. One of the main characteristics of the cat state is the superposition of two coherent states with opposite phase.
 
diff --git a/QuantumToolbox.jl/time_evolution/lowrank.qmd b/QuantumToolbox.jl/time_evolution/lowrank.qmd
index 32f5902..62a46a1 100644
--- a/QuantumToolbox.jl/time_evolution/lowrank.qmd
+++ b/QuantumToolbox.jl/time_evolution/lowrank.qmd
@@ -1,7 +1,7 @@
 ---
 title: Low rank master equation
 author: Luca Gravina
-date: 2025-05-22  # last update (keep this comment as a reminder)
+date: last-modified
 ---
 
 ## Introduction
@@ -33,6 +33,8 @@ In this example we consider the dynamics of the transverse field Ising model (TF
 using QuantumToolbox
 using LinearAlgebra
 using CairoMakie
+
+CairoMakie.activate!(type = "svg")
 ```
 
 We define the lattice with dimensions `Nx = 2` and `Ny = 3` and use the `Lattice` class to generate the lattice.
@@ -219,7 +221,7 @@ We can now compare the results of the low-rank evolution with the full evolution
 m_me = real(sol_me.expect[3, :]) / Nx / Ny
 m_lr = real(sol_lr.expect[3, :]) / Nx / Ny
 
-fig = Figure(size = (500, 350), fontsize = 15)
+fig = Figure(size = (900, 400), fontsize = 15)
 ax = Axis(fig[1, 1], xlabel = L"\gamma t", ylabel = L"M_{z}", xlabelsize = 20, ylabelsize = 20)
 lines!(ax, tl, m_lr, label = L"LR $[M=M(t)]$", linewidth = 2)
 lines!(ax, tl, m_me, label = "Fock", linewidth = 2, linestyle = :dash)
diff --git a/QuantumToolbox.jl/time_evolution/rabi.qmd b/QuantumToolbox.jl/time_evolution/rabi.qmd
index 40ae3d3..241c1e2 100644
--- a/QuantumToolbox.jl/time_evolution/rabi.qmd
+++ b/QuantumToolbox.jl/time_evolution/rabi.qmd
@@ -1,7 +1,7 @@
 ---
 title: "Vacuum Rabi oscillation"
 author: Li-Xun Cai
-date: 2025-05-22 # last update (keep this comment as a reminder)
+date: last-modified
 ---
     
 Inspirations taken from [this QuTiP tutorial](https://nbviewer.org/urls/qutip.org/qutip-tutorials/tutorials-v5/time-evolution/004_rabi-oscillations.ipynb) by J.R. Johansson, P.D. Nation, and C. Staufenbiel
@@ -42,9 +42,10 @@ $$
 
 ##### import:
 ```{julia}
-import QuantumToolbox: sigmaz, sigmam, destroy, qeye, basis, fock, 
-    ⊗, n_thermal, sesolve, mesolve
-import CairoMakie: Figure, Axis, lines!, axislegend, xlims!, display, @L_str
+using QuantumToolbox
+using CairoMakie
+
+CairoMakie.activate!(type = "svg")
 ```
 
 ```{julia}
@@ -84,28 +85,26 @@ eop_ls = [
 ]
 
 sol = sesolve(Htot , ψ0, tlist; e_ops = eop_ls)
-print(sol)
+sol
 ```
 
 Compare the dynamics of $| e \rangle\langle e|$ alongside $a^\dagger a$
 ```{julia}
 n = real.(sol.expect[1, :])
 e = real.(sol.expect[2, :])
-fig_se = Figure(size = (600, 350))
+fig_se = Figure()
 ax_se = Axis(
     fig_se[1, 1],
     xlabel = L"time $[1/\omega_a]$", 
     ylabel = "expectation value", 
     xlabelsize = 15, 
     ylabelsize = 15,
-    width = 400,
-    height = 220
 )
 xlims!(ax_se, 0, 400)
 lines!(ax_se, tlist, n, label = L"$\langle a^\dagger a \rangle$")
 lines!(ax_se, tlist, e, label = L"$P_e$")
 axislegend(ax_se; position = :rt, labelsize = 15)
-display(fig_se);
+fig_se
 ```
 
 In the above plot, the behaviour of the energy exchange between the atom and the cavity is clearly visible, addressing the Rabi problem.
@@ -180,20 +179,18 @@ print(sol_me)
 n_me = real.(sol_me.expect[1, :])
 e_me = real.(sol_me.expect[2, :])
 
-fig_me = Figure(size = (600, 350))
+fig_me = Figure()
 ax_me = Axis(
     fig_me[1, 1],
     xlabel = L"time $[1/\omega_a]$", 
     ylabel = "expectation value", 
     xlabelsize = 15, 
     ylabelsize = 15,
-    width = 400,
-    height = 220
 )
 lines!(ax_me, tlist, n_me, label = L"\langle a^\dagger a \rangle")
 lines!(ax_me, tlist, e_me, label = L"$P_e$")
 axislegend(ax_me; position = :rt, labelsize = 15)
-display(fig_me);
+fig_me
 ```
 
 From the above example, one can see that the dissipative system is losing energy over time and asymptoting to zero. We can further consider the thermal field with finite temperature.
@@ -203,20 +200,18 @@ sol_me_  = mesolve(Htot,  ψ0, tlist, cop_ls(γ, κ, 0.3 * ωa), e_ops = eop_ls)
 
 n_me_ = real.(sol_me_.expect[1, :])
 e_me_ = real.(sol_me_.expect[2, :])
-fig_me_ = Figure(size = (600, 350))
+fig_me_ = Figure()
 ax_me_ = Axis(
     fig_me_[1, 1],
     xlabel = L"time $[1/\omega_a]$", 
     ylabel = "expectation value", 
     xlabelsize = 15, 
     ylabelsize = 15,
-    width = 400,
-    height = 220
 )
 lines!(ax_me_, tlist, n_me_, label = L"\langle a^\dagger a \rangle")
 lines!(ax_me_, tlist, e_me_, label = L"$P_e$")
 axislegend(ax_me_; position = :rt, labelsize = 15)
-display(fig_me_);
+fig_me_
 ```
 Despite the persistence of the asymptotic behaviour, the system no longer approaches zero but instead reaches a steady-state above zero. This indicates that the system has been thermalized by the environment.
 
diff --git a/_quarto.yml b/_quarto.yml
index 2ce9184..60f296b 100644
--- a/_quarto.yml
+++ b/_quarto.yml
@@ -63,8 +63,10 @@ format:
     footnotes-hover: true
     fig-align: center
     fig-responsive: true
+    fig-width: 8.5
+    fig-height: 4
 
 execute:
   eval: true
-  cache: true
+  cache: false
   freeze: false