diff --git a/.buildkite/CUDA_Ext.yml b/.buildkite/CUDA_Ext.yml
index 40a80d6e5..6db9c208b 100644
--- a/.buildkite/CUDA_Ext.yml
+++ b/.buildkite/CUDA_Ext.yml
@@ -1,9 +1,41 @@
 steps:
+  # Julia 1.10 doesn't support [sources] in Project.toml,
+  # so we need to manually dev workspace packages before instantiating.
+  - label: "CUDA Julia 1.10"
+    plugins:
+      - JuliaCI/julia#v1:
+          version: "1.10"
+      - JuliaCI/julia-coverage#v1:
+          codecov: true
+          dirs:
+            - src
+            - lib
+            - ext
+    agents:
+      queue: "juliagpu"
+      cuda: "*"
+    commands: |
+      julia --project=. -e '
+        using Pkg
+        if VERSION < v"1.11"
+          for lib in readdir("lib"; join=true)
+            isdir(lib) || continue
+            isfile(joinpath(lib, "Project.toml")) || continue
+            Pkg.develop(path=lib)
+          end
+        end
+        Pkg.test(; test_args=["--quickfail"])
+      '
+    env:
+      GROUP: "CUDA-Ext"
+      SECRET_CODECOV_TOKEN: "ZfhQu/IcRLqNyZ//ZNs5sjBPaV76IHfU5gui52Qn+Rp8tOurukqgScuyDt+3HQ4R0hJYBw1/Nqg6jmBsvWSc9NEUx8kGsUJFHfN3no0+b+PFxA8oJkWc9EpyIsjht5ZIjlsFWR3f0DpPqMEle/QyWOPcal63CChXR8oAoR+Fz1Bh8GkokLlnC8F9Ugp9xBlu401GCbyZhvLTZnNIgK5yy9q8HBJnBg1cPOhI81J6JvYpEmcIofEzFV/qkfpTUPclu43WNoFX2DZPzbxilf3fsAd5/+nRkRfkNML8KiN4mnmjHxPPbuY8F5zC/PS5ybXtDpfvaMQc01WApXCkZk0ZAQ==;U2FsdGVkX1+eDT7dqCME5+Ox5i8GvWRTQbwiP/VYjapThDbxXFDeSSIC6Opmon+M8go22Bun3bat6Fzie65ang=="
+    timeout_in_minutes: 60
+
   - label: "CUDA Julia {{matrix.version}}"
     matrix:
       setup:
         version:
-          - "1.10" # oldest
+          - "1.11" # oldest
           - "1"    # latest
     plugins:
       - JuliaCI/julia#v1:
@@ -14,11 +46,12 @@ steps:
           codecov: true
           dirs:
             - src
+            - lib
             - ext
     agents:
       queue: "juliagpu"
       cuda: "*"
     env:
-      GROUP: "CUDA_Ext"
+      GROUP: "CUDA-Ext"
       SECRET_CODECOV_TOKEN: "ZfhQu/IcRLqNyZ//ZNs5sjBPaV76IHfU5gui52Qn+Rp8tOurukqgScuyDt+3HQ4R0hJYBw1/Nqg6jmBsvWSc9NEUx8kGsUJFHfN3no0+b+PFxA8oJkWc9EpyIsjht5ZIjlsFWR3f0DpPqMEle/QyWOPcal63CChXR8oAoR+Fz1Bh8GkokLlnC8F9Ugp9xBlu401GCbyZhvLTZnNIgK5yy9q8HBJnBg1cPOhI81J6JvYpEmcIofEzFV/qkfpTUPclu43WNoFX2DZPzbxilf3fsAd5/+nRkRfkNML8KiN4mnmjHxPPbuY8F5zC/PS5ybXtDpfvaMQc01WApXCkZk0ZAQ==;U2FsdGVkX1+eDT7dqCME5+Ox5i8GvWRTQbwiP/VYjapThDbxXFDeSSIC6Opmon+M8go22Bun3bat6Fzie65ang=="
     timeout_in_minutes: 60
diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml
index 36199fd77..b4e69a8dd 100644
--- a/.buildkite/pipeline.yml
+++ b/.buildkite/pipeline.yml
@@ -18,6 +18,7 @@ steps:
             - ".buildkite/pipeline.yml"
             - ".buildkite/CUDA_Ext.yml"
             - "src/**"
+            - "lib/**"
             - "ext/QuantumToolboxCUDAExt.jl"
             - "test/runtests.jl"
             - "test/ext-test/gpu/**"
diff --git a/.github/copilot-instructions.md b/.github/copilot-instructions.md
index b18214474..4a7028bb4 100644
--- a/.github/copilot-instructions.md
+++ b/.github/copilot-instructions.md
@@ -82,11 +82,11 @@ QuantumToolbox.settings.ProgressMeterKWARGS = (showspeed=true, printed=true)
 ### Running Tests
 ```bash
 make test                    # Run all core tests
-GROUP=Core make test         # Run core tests only
-GROUP=CUDA_Ext make test     # Run CUDA extension tests only
+GROUP=Main make test         # Run Main tests only
+GROUP=CUDA-Ext make test     # Run CUDA extension tests only
 ```
 
-Test groups (set `GROUP` env var): "All", "Core", "Code-Quality", "AutoDiff_Ext", "Makie_Ext", "CUDA_Ext", "Arbitrary-Precision"
+Test groups (set `GROUP` env var): "All", "Main", "Code-Quality", "AutoDiff-Ext", "Makie-Ext", "CUDA-Ext", "Arbitrary-Precision"
 
 ### Code Quality & Formatting
 ```bash
diff --git a/.github/workflows/Benchmarks.yml b/.github/workflows/Benchmarks.yml
index 5592d8c27..d73c8b337 100644
--- a/.github/workflows/Benchmarks.yml
+++ b/.github/workflows/Benchmarks.yml
@@ -6,6 +6,7 @@ on:
     paths:
       - '.github/workflows/Benchmarks.yml'
       - 'src/**'
+      - 'lib/**'
       - 'ext/**'
       - 'benchmarks/**'
       - 'Project.toml'
@@ -15,6 +16,7 @@ on:
     paths:
       - '.github/workflows/Benchmarks.yml'
       - 'src/**'
+      - 'lib/**'
       - 'ext/**'
       - 'benchmarks/**'
       - 'Project.toml'
diff --git a/.github/workflows/CI-Julia-nightly.yml b/.github/workflows/CI-Julia-nightly.yml
index dede3beb5..03c87188b 100644
--- a/.github/workflows/CI-Julia-nightly.yml
+++ b/.github/workflows/CI-Julia-nightly.yml
@@ -7,6 +7,7 @@ on:
     paths:
       - '.github/workflows/CI-Julia-nightly.yml'
       - 'src/**'
+      - 'lib/**'
       - 'ext/**'
       - 'test/runtests.jl'
       - 'test/core-test/**'
@@ -17,6 +18,7 @@ on:
     paths:
       - '.github/workflows/CI-Julia-nightly.yml'
       - 'src/**'
+      - 'lib/**'
       - 'ext/**'
       - 'test/runtests.jl'
       - 'test/core-test/**'
@@ -45,7 +47,7 @@ jobs:
         arch:
           - 'x64'
         group:
-          - 'Core'
+          - 'Main'
 
     steps:
       - uses: actions/checkout@v6
diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index 791b4e47f..9e840ea8c 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -7,6 +7,7 @@ on:
     paths:
       - '.github/workflows/CI.yml'
       - 'src/**'
+      - 'lib/**'
       - 'ext/**'
       - 'test/runtests.jl'
       - 'test/core-test/**'
@@ -18,6 +19,7 @@ on:
     paths:
       - '.github/workflows/CI.yml'
       - 'src/**'
+      - 'lib/**'
       - 'ext/**'
       - 'test/runtests.jl'
       - 'test/core-test/**'
@@ -52,7 +54,7 @@ jobs:
           - os: 'macOS-latest'
             arch: 'arm64'
         group:
-          - 'Core'
+          - 'Main'
 
         include:
           # for core tests (intermediate versions)
@@ -60,19 +62,19 @@ jobs:
             node:
               os: 'ubuntu-latest'
               arch: 'x64'
-            group: 'Core'
+            group: 'Main'
 
           # for extension tests
           - version: '1'
             node:
               os: 'ubuntu-latest'
               arch: 'x64'
-            group: 'Makie_Ext'
+            group: 'Makie-Ext'
           - version: '1.12'
             node:
               os: 'ubuntu-latest'
               arch: 'x64'
-            group: 'AutoDiff_Ext'
+            group: 'AutoDiff-Ext'
           - version: '1'
             node:
               os: 'ubuntu-latest'
@@ -86,6 +88,20 @@ jobs:
           version: ${{ matrix.version }}
           arch: ${{ matrix.node.arch }}
       - uses: julia-actions/cache@v3
+      - name: Dev local libs for Julia < 1.11
+        shell: bash
+        run: |
+          julia --project=. -e '
+            using Pkg
+            if VERSION < v"1.11"
+              for lib in readdir("lib"; join=true)
+                isdir(lib) || continue
+                isfile(joinpath(lib, "Project.toml")) || continue
+                Pkg.develop(path=lib)
+              end
+            end
+            Pkg.instantiate()
+          '
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
         env:
diff --git a/.github/workflows/Code-Quality.yml b/.github/workflows/Code-Quality.yml
index f6a4aba9a..d6a293144 100644
--- a/.github/workflows/Code-Quality.yml
+++ b/.github/workflows/Code-Quality.yml
@@ -7,6 +7,7 @@ on:
     paths:
       - '.github/workflows/Code-Quality.yml'
       - 'src/**'
+      - 'lib/**'
       - 'ext/**'
       - 'test/runtests.jl'
       - 'test/core-test/**'
@@ -17,9 +18,10 @@ on:
     paths:
       - '.github/workflows/Code-Quality.yml'
       - 'src/**'
+      - 'lib/**'
       - 'ext/**'
       - 'test/runtests.jl'
-      - 'test/core-test/**'
+      - 'test/code-quality/**'
       - 'Project.toml'
     types:
       - opened
@@ -39,14 +41,12 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - 'lts'
+          - '1.10' # oldest
           - '1'
         os:
           - 'ubuntu-latest'
         arch:
           - 'x64'
-        group:
-          - 'Code-Quality'
 
     steps:
       - uses: actions/checkout@v6
@@ -55,7 +55,21 @@ jobs:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
       - uses: julia-actions/cache@v3
+      - name: Dev local libs for Julia < 1.11
+        shell: bash
+        run: |
+          julia --project=. -e '
+            using Pkg
+            if VERSION < v"1.11"
+              for lib in readdir("lib"; join=true)
+                isdir(lib) || continue
+                isfile(joinpath(lib, "Project.toml")) || continue
+                Pkg.develop(path=lib)
+              end
+            end
+            Pkg.instantiate()
+          '
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
         env:
-          GROUP: ${{ matrix.group }}
+          GROUP: 'Code-Quality'
diff --git a/.github/workflows/bump-version.yml b/.github/workflows/bump-version.yml
index 9b292ba4e..e215effa1 100644
--- a/.github/workflows/bump-version.yml
+++ b/.github/workflows/bump-version.yml
@@ -3,6 +3,14 @@ name: Bump Version
 on:
   workflow_dispatch:
     inputs:
+      package:
+        description: 'Package to bump'
+        required: true
+        default: 'QuantumToolbox'
+        type: choice
+        options:
+        - QuantumToolbox
+        - QuantumToolboxUtils
       bump_type:
         description: 'Version bump type'
         required: true
@@ -37,12 +45,20 @@ jobs:
           julia -e '
             using Pkg; Pkg.add("Changelog")
             using TOML, Dates, Changelog
-            project_file = "Project.toml"
+            package = ENV["PACKAGE"]
+            project_file = package == "QuantumToolbox" ? "Project.toml" : "lib/$(package)/Project.toml"
             changelog_file = "CHANGELOG.md"
 
             # --- PART 1: Calculate new version number ---
+            # the current version of QuantumToolbox 
+            p_QT = TOML.parsefile("Project.toml")
+            v_QT = VersionNumber(p_QT["version"])
+            
+            # the current version of selected package 
             p = TOML.parsefile(project_file)
             v = VersionNumber(p["version"])
+
+            # the version bump type (major, minor, patch)
             type = ENV["BUMP_TYPE"]
             new_v = if type == "major"
                 VersionNumber(v.major + 1, 0, 0)
@@ -53,9 +69,10 @@ jobs:
             end
 
             # --- PART 2: Output version number to GitHub Actions variables ---
-            println("Bump version from v$v to v$new_v")
+            println("Bump $package version from v$v to v$new_v")
             open(ENV["GITHUB_OUTPUT"], "a") do io
               println(io, "new_version=$new_v")
+              println(io, "project_file=$project_file")
             end
 
             # --- PART 3: Update Project.toml (using Regex to preserve formatting) ---
@@ -66,35 +83,79 @@ jobs:
 
             # --- PART 4: Update CHANGELOG.md ---
             content = read(changelog_file, String)
-            marker = "## [Unreleased](https://github.com/qutip/QuantumToolbox.jl/tree/main)"
-            insertion = "$marker\n\n\n\n## [v$new_v]\nRelease date: $(Dates.format(now(), "yyyy-mm-dd"))"
-            new_content = replace(content, marker => insertion)
+
+            marker1 = "\n\n## [v$v_QT]"
+            insertion1 = "\n- [lib] Bump `$package` version to `v$new_v`.\n$marker1"
+            new_content = replace(content, marker1 => insertion1)
+
+            if package == "QuantumToolbox"
+              marker2 = "## [Unreleased](https://github.com/qutip/QuantumToolbox.jl/tree/main)"
+              insertion2 = "$marker2\n\n\n\n## [v$new_v]\nRelease date: $(Dates.format(now(), "yyyy-mm-dd"))"
+              new_content = replace(new_content, marker2 => insertion2)
+            end
+
             write(changelog_file, new_content)
             Changelog.generate(Changelog.CommonMark(), changelog_file; repo = "qutip/QuantumToolbox.jl")
             println("Updated $changelog_file")
           '
         env:
+          PACKAGE: ${{ inputs.package }}
           BUMP_TYPE: ${{ inputs.bump_type }}
 
       - name: Create Branch
         run: |
           git config user.name "github-actions[bot]"
           git config user.email "github-actions[bot]@users.noreply.github.com"
-          
-          BRANCH_NAME="bump-to-${{ steps.updates.outputs.new_version }}"
+
+          BRANCH_NAME="bump-${{ inputs.package }}-version-to-v${{ steps.updates.outputs.new_version }}"
           git checkout -b $BRANCH_NAME
-          git add Project.toml CHANGELOG.md
-          git commit -m "[no ci] Bump version to v${{ steps.updates.outputs.new_version }}"
+          git add "${{ steps.updates.outputs.project_file }}" CHANGELOG.md
+          git commit -m "[no ci] Bump ${{ inputs.package }}.jl version to v${{ steps.updates.outputs.new_version }}"
           git push origin $BRANCH_NAME
 
       - name: Create Pull Request
+        id: create_pr
         env:
           GH_TOKEN: ${{ secrets.GITHUB_TOKEN }}
         run: |
-          BRANCH_NAME="bump-to-${{ steps.updates.outputs.new_version }}"
-          gh pr create \
+          BRANCH_NAME="bump-${{ inputs.package }}-version-to-v${{ steps.updates.outputs.new_version }}"
+          PR_URL=$(gh pr create \
             --base main \
             --head $BRANCH_NAME \
-            --title "Bump version to v${{ steps.updates.outputs.new_version }}" \
-            --body "This PR bumps the package version to v${{ steps.updates.outputs.new_version }} and updates CHANGELOG.md." \
-            --label "Skip ChangeLog"
+            --title "Bump ${{ inputs.package }} version to v${{ steps.updates.outputs.new_version }}" \
+            --body "This PR bumps ${{ inputs.package }} version to v${{ steps.updates.outputs.new_version }} and updates CHANGELOG.md.")
+
+          PR_NUMBER="${PR_URL##*/}"
+
+          echo "Created PR: $PR_URL"
+          {
+            echo "pr_url=$PR_URL"
+            echo "pr_number=$PR_NUMBER"
+          } >> "$GITHUB_OUTPUT"
+
+      - name: Update CHANGELOG with PR number
+        run: |
+          julia -e '
+            using Pkg; Pkg.add("Changelog")
+            using Changelog
+
+            package = ENV["PACKAGE"]
+            new_v = ENV["NEW_VERSION"]
+            pr_number = ENV["PR_NUMBER"]
+            changelog_file = "CHANGELOG.md"
+
+            content = read(changelog_file, String)
+            old_entry = "- [lib] Bump `$package` version to `v$new_v`."
+            new_entry = "$old_entry ([#$pr_number])"
+
+            write(changelog_file, replace(content, old_entry => new_entry))
+            Changelog.generate(Changelog.CommonMark(), changelog_file; repo = "qutip/QuantumToolbox.jl")
+          '
+
+          git add CHANGELOG.md
+          git commit -m "[no ci] Add PR reference to CHANGELOG"
+          git push origin "bump-${{ inputs.package }}-version-to-v${{ steps.updates.outputs.new_version }}"
+        env:
+          PACKAGE: ${{ inputs.package }}
+          NEW_VERSION: ${{ steps.updates.outputs.new_version }}
+          PR_NUMBER: ${{ steps.create_pr.outputs.pr_number }}
diff --git a/.github/workflows/documentation.yml b/.github/workflows/documentation.yml
index 44b46e391..141240a6f 100644
--- a/.github/workflows/documentation.yml
+++ b/.github/workflows/documentation.yml
@@ -12,6 +12,7 @@ on:
     paths:
       - '.github/workflows/documentation.yml'
       - 'src/**'
+      - 'lib/**'
       - 'ext/**'
       - 'docs/**'
       - 'Project.toml'
diff --git a/CHANGELOG.md b/CHANGELOG.md
index de9f3d161..9c1cd656f 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -7,7 +7,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased](https://github.com/qutip/QuantumToolbox.jl/tree/main)
 
-
+- [lib] Introduce QuantumToolboxUtils library. ([#686])
 
 ## [v0.45.0]
 Release date: 2026-04-19
@@ -484,6 +484,7 @@ Release date: 2024-11-13
 [#681]: https://github.com/qutip/QuantumToolbox.jl/issues/681
 [#682]: https://github.com/qutip/QuantumToolbox.jl/issues/682
 [#683]: https://github.com/qutip/QuantumToolbox.jl/issues/683
+[#686]: https://github.com/qutip/QuantumToolbox.jl/issues/686
 [#688]: https://github.com/qutip/QuantumToolbox.jl/issues/688
 [#689]: https://github.com/qutip/QuantumToolbox.jl/issues/689
 [#690]: https://github.com/qutip/QuantumToolbox.jl/issues/690
diff --git a/Makefile b/Makefile
index 7e59b10a7..0ef7cc53e 100644
--- a/Makefile
+++ b/Makefile
@@ -13,11 +13,10 @@ changelog:
 	${JULIA} -e 'using Changelog; Changelog.generate(Changelog.CommonMark(), "CHANGELOG.md"; repo = "qutip/QuantumToolbox.jl")'
 
 test:
-	${JULIA} --project -e 'using Pkg; Pkg.update(); Pkg.test()'
+	${JULIA} -t 4 --project -e 'using Pkg; Pkg.update(); Pkg.test()'
 
 docs:
-	${JULIA} --project=docs -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.update()'
-	${JULIA} --project=docs docs/make.jl
+	${JULIA} --project=docs -e 'using Pkg; Pkg.update(); include("docs/make.jl")'
 
 vitepress:
 	npm --prefix docs run docs:dev
diff --git a/Project.toml b/Project.toml
index d77c82a89..31c715ecd 100644
--- a/Project.toml
+++ b/Project.toml
@@ -20,7 +20,9 @@ OrdinaryDiffEqLowOrderRK = "1344f307-1e59-4825-a18e-ace9aa3fa4c6"
 OrdinaryDiffEqVerner = "79d7bb75-1356-48c1-b8c0-6832512096c2"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
+QuantumToolboxUtils = "02afad9e-e2b8-4f9a-965f-b5dca5b538ca"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
@@ -36,6 +38,9 @@ GPUArrays = "0c68f7d7-f131-5f86-a1c3-88cf8149b2d7"
 KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
 Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 
+[sources]
+QuantumToolboxUtils = {path = "lib/QuantumToolboxUtils"}
+
 [extensions]
 QuantumToolboxCUDAExt = "CUDA"
 QuantumToolboxChainRulesCoreExt = "ChainRulesCore"
@@ -64,7 +69,9 @@ OrdinaryDiffEqLowOrderRK = "1"
 OrdinaryDiffEqVerner = "1"
 Pkg = "1"
 ProgressMeter = "1.11.0"
+QuantumToolboxUtils = "0.1"
 Random = "1"
+Reexport = "1"
 SciMLBase = "2.105"
 SciMLOperators = "1.12"
 SparseArrays = "1"
diff --git a/docs/Project.toml b/docs/Project.toml
index 57a8c133b..6af8c6a5c 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -8,7 +8,12 @@ DocumenterVitepress = "4710194d-e776-4893-9690-8d956a29c365"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Mooncake = "da2b9cff-9c12-43a0-ae48-6db2b0edb7d6"
 QuantumToolbox = "6c2fb7c5-b903-41d2-bc5e-5a7c320b9fab"
+QuantumToolboxUtils = "02afad9e-e2b8-4f9a-965f-b5dca5b538ca"
 SciMLSensitivity = "1ed8b502-d754-442c-8d5d-10ac956f44a1"
 
+[sources]
+QuantumToolbox = {path = ".."}
+QuantumToolboxUtils = {path = "../lib/QuantumToolboxUtils"}
+
 [compat]
 DocumenterVitepress = "0.3"
diff --git a/docs/make.jl b/docs/make.jl
index 59ab77ec9..0d0f07040 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -7,17 +7,15 @@ using Changelog
 # Load of packages required to compile the extension documentation
 using CairoMakie
 
-doctest_setup = quote
-    using LinearAlgebra
-    using SparseArrays
-    using QuantumToolbox
-end
-DocMeta.setdocmeta!(QuantumToolbox, :DocTestSetup, doctest_setup; recursive = true)
-
 # some options for `makedocs`
 const DRAFT = get(ENV, "DRAFT", false) == "true"  # `DRAFT   = true`  disables cell evaluation
 const DOCTEST = get(ENV, "DOCTEST", true) == true # `DOCTEST = false` skips doc tests
 
+if DOCTEST
+    DocMeta.setdocmeta!(QuantumToolboxUtils, :DocTestSetup, :(using QuantumToolboxUtils); recursive = true)
+    DocMeta.setdocmeta!(QuantumToolbox, :DocTestSetup, :(using QuantumToolbox, LinearAlgebra, SparseArrays); recursive = true)
+end
+
 # generate bibliography
 bib = CitationBibliography(
     joinpath(@__DIR__, "src", "resources", "bibliography.bib"),
@@ -36,7 +34,6 @@ const PAGES = [
     "Home" => "index.md",
     "Getting Started" => [
         "Brief Example" => "getting_started/brief_example.md",
-        # "Key differences from QuTiP" => "getting_started/qutip_differences.md",
         "The Importance of Type-Stability" => "getting_started/type_stability.md",
         "Example: Create QuantumToolbox.jl Logo" => "getting_started/logo.md",
         "Cite QuantumToolbox.jl" => "getting_started/cite.md",
@@ -83,6 +80,7 @@ const PAGES = [
 
 makedocs(;
     modules = [
+        QuantumToolboxUtils,
         QuantumToolbox,
         Base.get_extension(QuantumToolbox, :QuantumToolboxMakieExt),
     ],
diff --git a/docs/src/.vitepress/config.mts b/docs/src/.vitepress/config.mts
index 91a00ca79..1365173be 100644
--- a/docs/src/.vitepress/config.mts
+++ b/docs/src/.vitepress/config.mts
@@ -54,6 +54,10 @@ export default defineConfig({
       light: "github-light",
       dark: "github-dark"
     },
+
+    // options for @mdit-vue/plugin-toc
+    // https://github.com/mdit-vue/mdit-vue/tree/main/packages/plugin-toc#options
+    toc: { level: [2, 3, 4] }, // for API page, triggered by: [[toc]]
   },
   vite: {
     plugins: [
diff --git a/docs/src/index.md b/docs/src/index.md
index 154d27418..e0736c5a8 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -75,7 +75,7 @@ Alternatively, this can also be done in `Julia`'s [Pkg REPL](https://julialang.g
 ```
 More information about `Julia`'s package manager can be found at [`Pkg.jl`](https://julialang.github.io/Pkg.jl/v1/).
 
-To load the package and check the version information, use either [`QuantumToolbox.versioninfo()`](@ref) or [`QuantumToolbox.about()`](@ref), namely
+To load the package and check the version information, use either [`QuantumToolbox.versioninfo()`](@ref QuantumToolboxUtils.versioninfo) or [`QuantumToolbox.about()`](@ref QuantumToolboxUtils.about), namely
 ```julia
 using QuantumToolbox
 QuantumToolbox.versioninfo()
diff --git a/docs/src/resources/api.md b/docs/src/resources/api.md
index 435dbf91d..915d1660c 100644
--- a/docs/src/resources/api.md
+++ b/docs/src/resources/api.md
@@ -1,3 +1,34 @@
+# [API](@id doc-API)
+
+**Table of contents**
+
+[[toc]] <!-- the level setting is in ".vitepress/config.mts" -->
+
+## [QuantumToolboxUtils](@id doc-API:QuantumToolboxUtils)
+
+```@meta
+CurrentModule = QuantumToolboxUtils
+
+DocTestSetup = quote
+    using QuantumToolboxUtils
+end
+```
+
+```@docs
+settings
+versioninfo
+about
+cite
+meshgrid
+AbstractLinearMap
+gaussian
+n_thermal
+PhysicalConstants
+convert_unit
+```
+
+## [QuantumToolbox](@id doc-API:QuantumToolbox)
+
 ```@meta
 CurrentModule = QuantumToolbox
 
@@ -8,22 +39,23 @@ DocTestSetup = quote
 end
 ```
 
-# [API](@id doc-API)
-
-**Table of contents**
-
-[[toc]] <!-- the level setting is in ".vitepress/config.mts" -->
-
-## [Quantum object (Qobj) and type](@id doc-API:Quantum-object-and-type)
+### [Dimensions](@id doc-API:Dimensions)
 
 ```@docs
 Dimensions
+get_size
 AbstractSpace
 AbstractSuperSpace
 Space
 TensorSpace
 LiouvilleSpace
 EnrSpace
+enr_state_dictionaries
+```
+
+### [Quantum object (Qobj) and type](@id doc-API:Quantum-object-and-type)
+
+```@docs
 AbstractQuantumObject
 Bra
 Ket
@@ -39,7 +71,7 @@ Base.length
 SciMLOperators.cache_operator
 ```
 
-## [Qobj boolean functions](@id doc-API:Qobj-boolean-functions)
+### [Qobj boolean functions](@id doc-API:Qobj-boolean-functions)
 
 ```@docs
 isbra
@@ -57,7 +89,7 @@ SciMLOperators.iscached
 SciMLOperators.isconstant
 ```
 
-## [Qobj arithmetic and attributes](@id doc-API:Qobj-arithmetic-and-attributes)
+### [Qobj arithmetic and attributes](@id doc-API:Qobj-arithmetic-and-attributes)
 
 ```@docs
 Base.zero
@@ -89,7 +121,7 @@ get_coherence
 partial_transpose
 ```
 
-## [Qobj eigenvalues and eigenvectors](@id doc-API:Qobj-eigenvalues-and-eigenvectors)
+### [Qobj eigenvalues and eigenvectors](@id doc-API:Qobj-eigenvalues-and-eigenvectors)
 
 ```@docs
 EigsolveResult
@@ -101,7 +133,7 @@ eigsolve
 eigsolve_al
 ```
 
-## [Qobj manipulation](@id doc-API:Qobj-manipulation)
+### [Qobj manipulation](@id doc-API:Qobj-manipulation)
 
 ```@docs
 ket2dm
@@ -114,7 +146,7 @@ vec2mat
 mat2vec
 ```
 
-## [Generate states and operators](@id doc-API:Generate-states-and-operators)
+### [Generate states and operators](@id doc-API:Generate-states-and-operators)
 
 ```@docs
 zero_ket
@@ -172,7 +204,7 @@ sprepost
 lindblad_dissipator
 ```
 
-## [Synonyms of functions for Qobj](@id doc-API:Synonyms-of-functions-for-Qobj)
+### [Synonyms of functions for Qobj](@id doc-API:Synonyms-of-functions-for-Qobj)
 ```@docs
 Qobj
 QobjEvo
@@ -196,7 +228,7 @@ qeye_like
 qzero_like
 ```
 
-## [Time evolution](@id doc-API:Time-evolution)
+### [Time evolution](@id doc-API:Time-evolution)
 
 ```@docs
 TimeEvolutionProblem
@@ -229,7 +261,7 @@ brterm
 brmesolve
 ```
 
-### [Steady State Solvers](@id doc-API:Steady-State-Solvers)
+#### [Steady State Solvers](@id doc-API:Steady-State-Solvers)
 
 ```@docs
 steadystate
@@ -241,14 +273,14 @@ SteadyStateODESolver
 SSFloquetEffectiveLiouvillian
 ```
 
-### [Dynamical Shifted Fock method](@id doc-API:Dynamical-Shifted-Fock-method)
+#### [Dynamical Shifted Fock method](@id doc-API:Dynamical-Shifted-Fock-method)
 
 ```@docs
 dsf_mesolve
 dsf_mcsolve
 ```
 
-### [Low-rank time evolution](@id doc-API:Low-rank-time-evolution)
+#### [Low-rank time evolution](@id doc-API:Low-rank-time-evolution)
 
 ```@docs
 TimeEvolutionLRSol
@@ -256,7 +288,7 @@ lr_mesolveProblem
 lr_mesolve
 ```
 
-## [Correlations and Spectrum](@id doc-API:Correlations-and-Spectrum)
+### [Correlations and Spectrum](@id doc-API:Correlations-and-Spectrum)
 
 ```@docs
 correlation_3op_2t
@@ -270,7 +302,7 @@ PseudoInverse
 Lanczos
 ```
 
-## [Entropy and Metrics](@id doc-API:Entropy-and-Metrics)
+### [Entropy and Metrics](@id doc-API:Entropy-and-Metrics)
 
 ```@docs
 entropy_vn
@@ -289,7 +321,7 @@ bures_dist
 bures_angle
 ```
 
-## [Spin Lattice](@id doc-API:Spin-Lattice)
+### [Spin Lattice](@id doc-API:Spin-Lattice)
 
 ```@docs
 Lattice
@@ -297,45 +329,20 @@ multisite_operator
 DissipativeIsing
 ```
 
-## [Symmetries and Block Diagonalization](@id doc-API:Symmetries-and-Block-Diagonalization)
+### [Symmetries and Block Diagonalization](@id doc-API:Symmetries-and-Block-Diagonalization)
 
 ```@docs
 block_diagonal_form
 BlockDiagonalForm
 ```
 
-## [Miscellaneous](@id doc-API:Miscellaneous)
+### [Miscellaneous](@id doc-API:Miscellaneous)
 
 ```@docs
 wigner
 ```
 
-## [Linear Maps](@id doc-API:Linear-Maps)
-
-```@docs
-AbstractLinearMap
-```
-
-## [Utility functions](@id doc-API:Utility-functions)
-
-```@docs
-QuantumToolbox.settings
-QuantumToolbox.versioninfo
-QuantumToolbox.about
-QuantumToolbox.pkginfo
-QuantumToolbox.sysinfo
-QuantumToolbox.cite
-gaussian
-n_thermal
-PhysicalConstants
-convert_unit
-row_major_reshape
-get_size
-meshgrid
-enr_state_dictionaries
-```
-
-## [Visualization](@id doc-API:Visualization)
+### [Visualization](@id doc-API:Visualization)
 
 ```@docs
 plot_wigner
@@ -344,7 +351,7 @@ matrix_heatmap
 matrix_histogram
 ```
 
-### [Bloch Sphere](@id doc-API:Bloch-Sphere)
+#### [Bloch Sphere](@id doc-API:Bloch-Sphere)
 
 ```@docs
 Bloch
diff --git a/docs/src/resources/contributing.md b/docs/src/resources/contributing.md
index 494a19523..b255c3409 100644
--- a/docs/src/resources/contributing.md
+++ b/docs/src/resources/contributing.md
@@ -62,15 +62,15 @@ make test
 
 This command will automatically rebuild `Julia` and run the script located in `test/runtests.jl` (should cover both the original tests and the new test(s) you add).
 
-The tests are divided into several test groups, where the group names are defined in the file `test/runtests.jl` with a variable `GROUP`. One can also run the test scripts just for a certain test group by adding an argument `GROUP=<test-group-name>` to the `make test` command. For example, to run the tests for group `Core`, one can use the following command:
+The tests are divided into several test groups, where the group names are defined in the file `test/runtests.jl` with a variable `GROUP`. One can also run the test scripts just for a certain test group by adding an argument `GROUP=<test-group-name>` to the `make test` command. For example, to run the tests for group `Main`, one can use the following command:
 
 ```shell
-make GROUP=Core test
+make test GROUP=Main
 ```
 
-#### [Test Item Framework for Core tests](@id doc-Contribute:Test-Item-Framework-for-Core-tests)
+#### [Test Item Framework for Main tests](@id doc-Contribute:Test-Item-Framework-for-Main-tests)
 
-The tests in `GROUP=Core` are provided using the [Test Item Framework](https://www.julia-vscode.org/docs/stable/userguide/testitems/), which structures the test codes into `@testitems` and makes it easier to run individually.
+The tests in `GROUP=Main` are provided using the [Test Item Framework](https://www.julia-vscode.org/docs/stable/userguide/testitems/), which structures the test codes into `@testitems` and makes it easier to run individually.
 
 The [VS Code](https://code.visualstudio.com/) and its [Julia extension](https://www.julia-vscode.org/) provides us with options to run individual `@testitems`. It is much easier to find the specific core test that failed since the [Julia extension](https://www.julia-vscode.org/) in [VS Code](https://code.visualstudio.com/) will collect all core test failures and then display them in a structured way, directly at the place in the code where a specific core test failed. See [here](https://www.julia-vscode.org/docs/stable/userguide/testitems/) for more details.
 
diff --git a/docs/src/users_guide/settings.md b/docs/src/users_guide/settings.md
index 53e9f59c9..78b01a25e 100644
--- a/docs/src/users_guide/settings.md
+++ b/docs/src/users_guide/settings.md
@@ -2,7 +2,7 @@
 
 In this section, we introduce the default global settings used throughout the package and show how to modify them.
 
-All settings are stored in [`QuantumToolbox.settings`](@ref).
+All settings are stored in [`QuantumToolbox.settings`](@ref QuantumToolboxUtils.settings).
 
 !!! warning "Differences from QuTiP"
     Due to the differences in programming languages, solving algorithms, and many other reasons, these global settings (including their default values and usage) may be very different from those in `Python QuTiP`.
@@ -20,7 +20,7 @@ Here, we list out each setting along with the specific functions that will use i
 
 ## Change default settings
 
-First, we can check the current [`QuantumToolbox.settings`](@ref):
+First, we can check the current [`QuantumToolbox.settings`](@ref QuantumToolboxUtils.settings):
 
 ```@example settings
 using QuantumToolbox
diff --git a/lib/QuantumToolboxUtils/LICENSE b/lib/QuantumToolboxUtils/LICENSE
new file mode 100644
index 000000000..a3dbb37b8
--- /dev/null
+++ b/lib/QuantumToolboxUtils/LICENSE
@@ -0,0 +1,29 @@
+BSD 3-Clause License
+
+Copyright (c) 2022 to 2026 inclusive, QuTiP developers and contributors.
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without modification,
+are permitted provided that the following conditions are met:
+
+1. Redistributions of source code must retain the above copyright notice, this
+   list of conditions and the following disclaimer.
+
+2. Redistributions in binary form must reproduce the above copyright notice,
+   this list of conditions and the following disclaimer in the documentation
+   and/or other materials provided with the distribution.
+
+3. Neither the name of the copyright holder nor the names of its contributors
+   may be used to endorse or promote products derived from this software without
+   specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR
+ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
+ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
\ No newline at end of file
diff --git a/lib/QuantumToolboxUtils/Project.toml b/lib/QuantumToolboxUtils/Project.toml
new file mode 100644
index 000000000..69b9116e7
--- /dev/null
+++ b/lib/QuantumToolboxUtils/Project.toml
@@ -0,0 +1,19 @@
+name = "QuantumToolboxUtils"
+uuid = "02afad9e-e2b8-4f9a-965f-b5dca5b538ca"
+authors = ["Alberto Mercurio", "Yi-Te Huang"]
+version = "0.1.0"
+
+[deps]
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
+
+[compat]
+LinearAlgebra = "1"
+Random = "1"
+SciMLOperators = "1.12"
+SparseArrays = "1"
+StaticArraysCore = "1"
+julia = "1.10"
diff --git a/lib/QuantumToolboxUtils/src/PhysicalConstants.jl b/lib/QuantumToolboxUtils/src/PhysicalConstants.jl
new file mode 100644
index 000000000..f6730e229
--- /dev/null
+++ b/lib/QuantumToolboxUtils/src/PhysicalConstants.jl
@@ -0,0 +1,80 @@
+export PhysicalConstants, convert_unit
+
+@doc raw"""
+    const PhysicalConstants
+
+A `NamedTuple` which stores some constant values listed in [*CODATA recommended values of the fundamental physical constants: 2022*](https://physics.nist.gov/cuu/pdf/wall_2022.pdf)
+
+The current stored constants are:
+- `c` : (exact) speed of light in vacuum with unit ``[\textrm{m}\cdot\textrm{s}^{-1}]``
+- `G` : Newtonian constant of gravitation with unit ``[\textrm{m}^3\cdot\textrm{kg}^{−1}\cdot\textrm{s}^{−2}]``
+- `h` : (exact) Planck constant with unit ``[\textrm{J}\cdot\textrm{s}]``
+- `ħ` : reduced Planck constant with unit ``[\textrm{J}\cdot\textrm{s}]``
+- `e` : (exact) elementary charge with unit ``[\textrm{C}]``
+- `μ0` : vacuum magnetic permeability with unit ``[\textrm{N}\cdot\textrm{A}^{-2}]``
+- `ϵ0` : vacuum electric permittivity with unit ``[\textrm{F}\cdot\textrm{m}^{-1}]``
+- `k` : (exact) Boltzmann constant with unit ``[\textrm{J}\cdot\textrm{K}^{-1}]``
+- `NA` : (exact) Avogadro constant with unit ``[\textrm{mol}^{-1}]``
+
+# Examples
+
+```jldoctest
+julia> PhysicalConstants.ħ
+1.0545718176461565e-34
+```
+"""
+const PhysicalConstants = (
+    c = 299792458.0,
+    G = 6.6743e-11,
+    h = 6.62607015e-34,
+    ħ = 6.62607015e-34 / (2 * π),
+    e = 1.602176634e-19,
+    μ0 = 1.25663706127e-6,
+    ϵ0 = 8.8541878188e-12,
+    k = 1.380649e-23,
+    NA = 6.02214076e23,
+)
+
+# common energy units (the values below are all in the unit of Joule)
+const _energy_units::Dict{Symbol, Float64} = Dict(
+    :J => 1.0,
+    :eV => PhysicalConstants.e,
+    :meV => 1.0e-3 * PhysicalConstants.e,
+    :MHz => 1.0e6 * PhysicalConstants.h,
+    :GHz => 1.0e9 * PhysicalConstants.h,
+    :K => PhysicalConstants.k,
+    :mK => 1.0e-3 * PhysicalConstants.k,
+)
+
+@doc raw"""
+    convert_unit(value::Real, unit1::Symbol, unit2::Symbol)
+
+Convert the energy `value` from `unit1` to `unit2`. The `unit1` and `unit2` can be either the following `Symbol`:
+- `:J` : Joule
+- `:eV` : electron volt
+- `:meV` : milli-electron volt
+- `:MHz` : Mega-Hertz multiplied by Planck constant ``h``
+- `:GHz` : Giga-Hertz multiplied by Planck constant ``h``
+- `:K` : Kelvin multiplied by Boltzmann constant ``k``
+- `:mK` : milli-Kelvin multiplied by Boltzmann constant ``k``
+
+Note that we use the values stored in [`PhysicalConstants`](@ref) to do the conversion.
+
+# Examples
+
+```jldoctest
+julia> convert_unit(1, :eV, :J)
+1.602176634e-19
+
+julia> convert_unit(1, :GHz, :J)
+6.62607015e-25
+
+julia> round(convert_unit(1, :meV, :mK), digits=4)
+11604.5181
+```
+"""
+function convert_unit(value::T, unit1::Symbol, unit2::Symbol) where {T <: Real}
+    !haskey(_energy_units, unit1) && throw(ArgumentError("Invalid unit :$(unit1)"))
+    !haskey(_energy_units, unit2) && throw(ArgumentError("Invalid unit :$(unit2)"))
+    return _float_type(T)(value * (_energy_units[unit1] / _energy_units[unit2]))
+end
diff --git a/lib/QuantumToolboxUtils/src/QuantumToolboxUtils.jl b/lib/QuantumToolboxUtils/src/QuantumToolboxUtils.jl
new file mode 100644
index 000000000..45de9f7dc
--- /dev/null
+++ b/lib/QuantumToolboxUtils/src/QuantumToolboxUtils.jl
@@ -0,0 +1,32 @@
+module QuantumToolboxUtils
+
+using LinearAlgebra
+using SparseArrays
+import StaticArraysCore: SVector, MVector
+import Random: AbstractRNG
+
+import SciMLOperators: SciMLOperators, AbstractSciMLOperator
+
+## Basic utilities for QuantumToolbox libraries
+include("settings.jl")
+include("versioninfo.jl")
+include("type_handle.jl")
+
+## Linear Algebra
+include("linalg/linalg.jl")
+include("linalg/arnoldi.jl")
+include("linalg/linear_maps.jl")
+
+## Other physical related utilities
+include("PhysicalConstants.jl")
+include("physics_func.jl")
+
+## deprecated functions
+include("deprecated.jl")
+
+function __init__()
+    _register_qt_library!(QuantumToolboxUtils)
+    return nothing
+end
+
+end
diff --git a/lib/QuantumToolboxUtils/src/deprecated.jl b/lib/QuantumToolboxUtils/src/deprecated.jl
new file mode 100644
index 000000000..13af5ae49
--- /dev/null
+++ b/lib/QuantumToolboxUtils/src/deprecated.jl
@@ -0,0 +1,23 @@
+#=
+This file gathers all the deprecated names (structures, functions, or variables) which will be removed in the future major release.
+
+- Before the major release, the deprecated names will show warnings or just throw errors when they are called.
+- If the deprecated names were once exported, we will still export them here until next major release.
+- If we decide to push a major release, cleanup this file.
+
+Example 1 [throw errors if the deprecation is fundamental, and there is no replacement for it]:
+```
+export deprecated_foo
+deprecated_foo(args...; kwargs...) = error("`deprecated_foo` is deprecated and will be removed in next major release.")
+```
+
+Example 2 ["force" show warning and tell the users that there is a replacement for the deprecated function]:
+
+```
+export deprecated_foo
+function deprecated_foo(args...; kwargs...)
+    Base.depwarn("`deprecated_foo` is deprecated and will be removed in next major release, use `new_foo` instead.", :deprecated_foo, force = true)
+    new_foo(args...; kwargs...)
+end
+```
+=#
diff --git a/src/arnoldi.jl b/lib/QuantumToolboxUtils/src/linalg/arnoldi.jl
similarity index 100%
rename from src/arnoldi.jl
rename to lib/QuantumToolboxUtils/src/linalg/arnoldi.jl
diff --git a/lib/QuantumToolboxUtils/src/linalg/linalg.jl b/lib/QuantumToolboxUtils/src/linalg/linalg.jl
new file mode 100644
index 000000000..320fd72eb
--- /dev/null
+++ b/lib/QuantumToolboxUtils/src/linalg/linalg.jl
@@ -0,0 +1,12 @@
+export meshgrid
+
+@doc raw"""
+    meshgrid(x::AbstractVector, y::AbstractVector)
+
+Equivalent to [numpy meshgrid](https://numpy.org/doc/stable/reference/generated/numpy.meshgrid.html).
+"""
+function meshgrid(x::AbstractVector{T}, y::AbstractVector{T}) where {T}
+    X = reshape(repeat(x, inner = length(y)), length(y), length(x))
+    Y = repeat(y, outer = (1, length(x)))
+    return X, Y
+end
diff --git a/src/linear_maps.jl b/lib/QuantumToolboxUtils/src/linalg/linear_maps.jl
similarity index 90%
rename from src/linear_maps.jl
rename to lib/QuantumToolboxUtils/src/linalg/linear_maps.jl
index 93f550788..e5184165d 100644
--- a/src/linear_maps.jl
+++ b/lib/QuantumToolboxUtils/src/linalg/linear_maps.jl
@@ -28,7 +28,7 @@ It is typically represented as a matrix with dimensions given by `size`, and thi
 
 ## Example
 
-As an example, we now define the linear map used in the [`eigsolve_al`](@ref) function for Arnoldi-Lindblad eigenvalue solver:
+As an example, we now define the linear map used in the [`eigsolve_al`](@ref QuantumToolbox.eigsolve_al) function for Arnoldi-Lindblad eigenvalue solver:
 
 ```julia
 struct ArnoldiLindbladIntegratorMap{T,TS,TI} <: AbstractLinearMap{T,TS}
@@ -44,7 +44,7 @@ function LinearAlgebra.mul!(y::AbstractVector, A::ArnoldiLindbladIntegratorMap,
 end
 ```
 
-where `integrator` is the ODE integrator for the time-evolution. In this way, we can diagonalize this linear map using the [`eigsolve`](@ref) function.
+where `integrator` is the ODE integrator for the time-evolution. In this way, we can diagonalize this linear map using the [`eigsolve`](@ref QuantumToolbox.eigsolve) function.
 """
 abstract type AbstractLinearMap{T, TS} end
 
diff --git a/lib/QuantumToolboxUtils/src/physics_func.jl b/lib/QuantumToolboxUtils/src/physics_func.jl
new file mode 100644
index 000000000..7a37736de
--- /dev/null
+++ b/lib/QuantumToolboxUtils/src/physics_func.jl
@@ -0,0 +1,28 @@
+export gaussian, n_thermal
+
+@doc raw"""
+    gaussian(x::Number, μ::Number, σ::Number)
+
+Returns the gaussian function ``\exp \left[- \frac{(x - \mu)^2}{2 \sigma^2} \right]``,
+where ``\mu`` and ``\sigma^2`` are the mean and the variance respectively.
+"""
+gaussian(x::Number, μ::Number, σ::Number) = exp(-(x - μ)^2 / (2 * σ^2))
+
+@doc raw"""
+    n_thermal(ω::Real, ω_th::Real)
+
+Return the number of photons in thermal equilibrium for an harmonic oscillator mode with frequency ``\omega``, at the temperature described by ``\omega_{\textrm{th}} \equiv k_B T / \hbar``:
+```math
+n(\omega, \omega_{\textrm{th}}) = \frac{1}{e^{\omega/\omega_{\textrm{th}}} - 1},
+```
+where ``\hbar`` is the reduced Planck constant, and ``k_B`` is the Boltzmann constant.
+"""
+function n_thermal(ω::T1, ω_th::T2) where {T1 <: Real, T2 <: Real}
+    x = exp(ω / ω_th)
+    n = ((x != 1) && (ω_th > 0)) ? 1 / (x - 1) : 0
+    return _float_type(promote_type(T1, T2))(n)
+end
+
+_Ginibre_ensemble(rng::AbstractRNG, ::Type{T}, n::Int, rank::Int = n) where {T <: Complex} = randn(rng, T, n, rank) / sqrt(T(n))
+
+_Boltzmann_weight(β::T, E::Int) where {T <: Real} = (E != 0 || isfinite(β)) ? exp(-β * E) : one(T)
diff --git a/src/settings.jl b/lib/QuantumToolboxUtils/src/settings.jl
similarity index 83%
rename from src/settings.jl
rename to lib/QuantumToolboxUtils/src/settings.jl
index e4daa9c85..e439d326b 100644
--- a/src/settings.jl
+++ b/lib/QuantumToolboxUtils/src/settings.jl
@@ -17,12 +17,12 @@ end
 @doc raw"""
     QuantumToolbox.settings
 
-Contains all the default global settings of QuantumToolbox.jl.
+Contains all the default global settings of QuantumToolbox.jl and all related libraries.
 
 # List of settings
 
-- `tidyup_tol::Float64 = 1.0e-14` : tolerance for [`tidyup`](@ref), [`tidyup!`](@ref) and [`to_sparse`](@ref).
-- `auto_tidyup::Bool = true` : Automatically tidyup.
+- `tidyup_tol::Float64 = 1.0e-14` : tolerance for [`tidyup`](@ref QuantumToolbox.tidyup), [`tidyup!`](@ref QuantumToolbox.tidyup!) and [`to_sparse`](@ref QuantumToolbox.to_sparse).
+- `auto_tidyup::Bool = true` : Automatically `tidyup`.
 - `ProgressMeterKWARGS::NamedTuple = (showspeed = true, printed = true)` : Default keyword arguments for progress bar in [`ProgressMeter.jl`](https://github.com/timholy/ProgressMeter.jl). This allows the customization of progress bar.
 
 For detailed explanation of each settings, see our documentation [here](https://qutip.org/QuantumToolbox.jl/stable/users_guide/settings).
diff --git a/lib/QuantumToolboxUtils/src/type_handle.jl b/lib/QuantumToolboxUtils/src/type_handle.jl
new file mode 100644
index 000000000..8737c4034
--- /dev/null
+++ b/lib/QuantumToolboxUtils/src/type_handle.jl
@@ -0,0 +1,70 @@
+makeVal(x::Val{T}) where {T} = x
+makeVal(x) = Val(x)
+
+getVal(x::Val{T}) where {T} = T
+getVal(x) = x # getVal for any other type
+
+_non_static_array_warning(argname, arg::Tuple{}) =
+    throw(ArgumentError("The argument $argname must be a Tuple or a StaticVector of non-zero length."))
+_non_static_array_warning(argname, arg::Union{SVector{N, T}, MVector{N, T}, NTuple{N, T}}) where {N, T} = nothing
+_non_static_array_warning(argname, arg::AbstractVector{T}) where {T} =
+    @warn "The argument $argname should be a Tuple or a StaticVector for better performance. Try to use `$argname = $(Tuple(arg))` instead of `$argname = $arg`. " *
+    "Alternatively, you can do `import QuantumToolbox: SVector` " *
+    "and use `$argname = SVector(" *
+    join(arg, ", ") *
+    ")`." maxlog = 1
+
+# lazy tensor warning
+for AType in (:AbstractArray, :AbstractSciMLOperator)
+    for BType in (:AbstractArray, :AbstractSciMLOperator)
+        if AType == BType == :AbstractArray
+            @eval begin
+                _lazy_tensor_warning(::$AType, ::$BType) = nothing
+            end
+        else
+            @eval begin
+                _lazy_tensor_warning(A::$AType, B::$BType) =
+                    @warn "using lazy tensor (which can hurt performance) between data types: $(get_typename_wrapper(A)) and $(get_typename_wrapper(B))"
+            end
+        end
+    end
+end
+
+get_typename_wrapper(A) = Base.typename(typeof(A)).wrapper
+
+_dense_similar(A::AbstractArray, args...) = similar(A, args...)
+_dense_similar(A::AbstractSparseMatrix, args...) = similar(nonzeros(A), args...)
+
+_sparse_similar(A::AbstractArray, args...) = sparse(args...)
+
+# alias of abstract types
+const FloatOrComplex = Union{T, Complex{T}} where {T <: AbstractFloat}
+
+# functions for getting Float or Complex element type
+_float_type(::AbstractArray{T}) where {T <: Number} = _float_type(T)
+_float_type(::AbstractSciMLOperator{T}) where {T <: Number} = _float_type(T)
+_float_type(::Type{Int32}) = Float32
+_float_type(::Type{Int64}) = Float64
+_float_type(::Type{Complex{Int32}}) = Float32
+_float_type(::Type{Complex{Int64}}) = Float64
+_float_type(::Type{Complex{T}}) where {T <: Real} = T
+_float_type(T::Type{<:AbstractFloat}) = T # Allow other untracked Real types, like ForwardDiff.Dual
+_complex_float_type(::AbstractArray{T}) where {T <: Number} = _complex_float_type(T)
+_complex_float_type(::AbstractSciMLOperator{T}) where {T <: Number} = _complex_float_type(T)
+_complex_float_type(::Type{Int32}) = ComplexF32
+_complex_float_type(::Type{Int64}) = ComplexF64
+_complex_float_type(::Type{Float32}) = ComplexF32
+_complex_float_type(::Type{Float64}) = ComplexF64
+_complex_float_type(::Type{Complex{Int32}}) = ComplexF32
+_complex_float_type(::Type{Complex{Int64}}) = ComplexF64
+_complex_float_type(::Type{Complex{Float32}}) = ComplexF32
+_complex_float_type(::Type{Complex{Float64}}) = ComplexF64
+_complex_float_type(T::Type{<:AbstractFloat}) = Complex{T} # Allow other untracked Complex types, like ForwardDiff.Dual
+_complex_float_type(T::Type{<:Complex}) = T       # Allow other untracked Complex types, like ForwardDiff.Dual
+
+_convert_eltype_wordsize(::Type{T}, ::Val{64}) where {T <: Int} = Int64
+_convert_eltype_wordsize(::Type{T}, ::Val{32}) where {T <: Int} = Int32
+_convert_eltype_wordsize(::Type{T}, ::Val{64}) where {T <: AbstractFloat} = Float64
+_convert_eltype_wordsize(::Type{T}, ::Val{32}) where {T <: AbstractFloat} = Float32
+_convert_eltype_wordsize(::Type{Complex{T}}, ::Val{64}) where {T <: Union{Int, AbstractFloat}} = ComplexF64
+_convert_eltype_wordsize(::Type{Complex{T}}, ::Val{32}) where {T <: Union{Int, AbstractFloat}} = ComplexF32
diff --git a/src/versioninfo.jl b/lib/QuantumToolboxUtils/src/versioninfo.jl
similarity index 58%
rename from src/versioninfo.jl
rename to lib/QuantumToolboxUtils/src/versioninfo.jl
index 36a5fc584..82388420e 100644
--- a/src/versioninfo.jl
+++ b/lib/QuantumToolboxUtils/src/versioninfo.jl
@@ -1,77 +1,72 @@
 #=
-Command line output of information on QuantumToolbox, dependencies, and system information
+Reusable version information helpers for QuantumToolbox libraries.
 =#
 
-@doc raw"""
-    QuantumToolbox.versioninfo(io::IO=stdout)
+# Registry of all loaded QuantumToolbox libraries, populated via __init__ in each library.
+const _QT_LIBRARIES = Module[]
 
-Command line output of information on QuantumToolbox, dependencies, and system information, same as [`QuantumToolbox.about`](@ref).
+raw"""
+    _register_qt_library!(m::Module)
+
+Register a QuantumToolbox library module into the global registry. Each library should call this in its `__init__` function.
 """
-function versioninfo(io::IO = stdout)
-    # print introduction
-    println(
-        io,
-        "\n",
-        " QuantumToolbox.jl: Quantum Toolbox in Julia\n",
-        "≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡\n",
-        "Copyright © QuTiP team 2022 and later.\n",
-        "Current admin team:\n",
-        "    Alberto Mercurio and Yi-Te Huang\n",
-    )
+function _register_qt_library!(m::Module)
+    (m ∉ _QT_LIBRARIES) && pushfirst!(_QT_LIBRARIES, m) # use pushfirst! so that main API libraries are at the front of the registry (for better display order in versioninfo)
+    return nothing
+end
 
-    # print package information
-    pkg_tuple = (
-        QuantumToolbox,
-        SciMLOperators,
-        LinearSolve,
-        OrdinaryDiffEqCore,
-    )
-    pkginfo(io; pkgs = pkg_tuple)
+raw"""
+    _add_library_deps!(lib::Val, DEPpkgs::Vector{Module})
 
-    # print System information
-    sysinfo(io)
+Add new dependencies to `DEPpkgs` for the specified library. Extend this function in each QuantumToolbox library to declare its external (non-QuantumToolbox) dependencies.
+"""
+_add_library_deps!(lib::Val{T}, DEPpkgs::Vector{Module}) where {T} = throw(ArgumentError("Unknown QuantumToolbox library : $T"))
+function _add_library_deps!(lib::Val{:QuantumToolboxUtils}, DEPpkgs::Vector{Module})
+    _add_pkgs!(DEPpkgs, Module[SciMLOperators])
+    return nothing
+end
 
-    # print citation information
-    println(
-        io,
-        "+---------------------------------------------------+\n",
-        "| Please cite QuantumToolbox.jl in your publication |\n",
-        "+---------------------------------------------------+\n",
-        "For your convenience, a bibtex reference can be easily generated using `QuantumToolbox.cite()`.\n",
-    )
+# add pkgs2 into pkgs1 but ensure the uniqueness
+function _add_pkgs!(pkgs1::Vector{Module}, pkgs2::Vector{Module})
+    for pkg in unique(pkgs2)
+        (pkg ∉ pkgs1) && push!(pkgs1, pkg)
+    end
     return nothing
 end
 
-@doc raw"""
-    QuantumToolbox.pkginfo(io::IO=stdout; pkgs::NTuple{N, Module} = (QuantumToolbox,))
+raw"""
+    QuantumToolboxUtils.pkginfo(io::IO=stdout; pkgs::Vector{Module} = Module[])
 
-Command line output of version numbers for given tuple of packages: `pkgs`. Default to `(QuantumToolbox,)`.
+Command line output of version numbers for given vector of packages: `pkgs`.
 """
-function pkginfo(io::IO = stdout; pkgs::NTuple{N, Module} = (QuantumToolbox,)) where {N}
+function pkginfo(io::IO = stdout; pkgs::Vector{Module} = Module[], split_after::Union{Nothing, Int} = nothing)
     pkg_ver_list = map(pkgversion, pkgs)
     maxLen = max(5, maximum(length ∘ string, pkgs)) # maximum string length of package names (5 refer to "Julia")
 
+    separation_line = "------------------------------------"
     print(
         io,
         "Package information:\n",
         "====================================\n",
     )
     println(io, rpad("Julia", maxLen, " "), " Ver. ", VERSION) # print Julia version first
-    for (pkg, pkg_ver) in zip(pkgs, pkg_ver_list)
+    println(io, separation_line)
+    for (idx, (pkg, pkg_ver)) in enumerate(zip(pkgs, pkg_ver_list))
         println(io, rpad(pkg, maxLen, " "), " Ver. ", pkg_ver)
+        !isnothing(split_after) && (idx == split_after) && (idx < length(pkgs)) && println(io, separation_line)
     end
     print(io, "\n")
     return nothing
 end
 
-@doc raw"""
-    QuantumToolbox.sysinfo(io::IO=stdout)
+raw"""
+    QuantumToolboxUtils.sysinfo(io::IO=stdout)
 
 Command line output of system information.
 """
 function sysinfo(io::IO = stdout)
     cpu = Sys.cpu_info()
-    BLAS_info = BLAS.get_config().loaded_libs[1]
+    BLAS_info = LinearAlgebra.BLAS.get_config().loaded_libs[1]
     Sys.iswindows() ? OS_name = "Windows" : Sys.isapple() ? OS_name = "macOS" : OS_name = Sys.KERNEL
 
     println(
@@ -90,10 +85,50 @@ function sysinfo(io::IO = stdout)
     return nothing
 end
 
+function _print_versioninfo(io::IO = stdout)
+    println(
+        io,
+        "\n",
+        " QuantumToolbox.jl: Quantum Toolbox in Julia\n",
+        "≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡\n",
+        "Copyright © QuTiP team 2022 and later.\n",
+        "Current admin team:\n",
+        "    Alberto Mercurio and Yi-Te Huang\n",
+    )
+
+    DEPpkgs = _gen_dep_pkg_list()
+    pkginfo(io; pkgs = vcat(_QT_LIBRARIES, DEPpkgs), split_after = length(_QT_LIBRARIES))
+    sysinfo(io)
+
+    println(
+        io,
+        "+---------------------------------------------------+\n",
+        "| Please cite QuantumToolbox.jl in your publication |\n",
+        "+---------------------------------------------------+\n",
+        "For your convenience, a bibtex reference can be easily generated using `QuantumToolbox.cite()`.\n",
+    )
+    return nothing
+end
+
+function _gen_dep_pkg_list()
+    DEPpkgs = Module[]
+    for lib in _QT_LIBRARIES
+        _add_library_deps!(Val(nameof(lib)), DEPpkgs)
+    end
+    return DEPpkgs
+end
+
+@doc raw"""
+    QuantumToolbox.versioninfo(io::IO=stdout)
+
+Command line output of information on QuantumToolbox, dependencies, and system information, same as [`QuantumToolbox.about`](@ref QuantumToolboxUtils.about).
+"""
+versioninfo(io::IO = stdout) = _print_versioninfo(io)
+
 @doc raw"""
     QuantumToolbox.about(io::IO=stdout)
 
-Command line output of information on QuantumToolbox, dependencies, and system information, same as [`QuantumToolbox.versioninfo`](@ref).
+Command line output of information on QuantumToolbox, dependencies, and system information, same as [`QuantumToolbox.versioninfo`](@ref QuantumToolboxUtils.versioninfo).
 """
 about(io::IO = stdout) = versioninfo(io)
 
diff --git a/lib/QuantumToolboxUtils/test/PhysicalConstants.jl b/lib/QuantumToolboxUtils/test/PhysicalConstants.jl
new file mode 100644
index 000000000..eba203b55
--- /dev/null
+++ b/lib/QuantumToolboxUtils/test/PhysicalConstants.jl
@@ -0,0 +1,38 @@
+@testitem "CODATA Physical Constants" begin
+    c = PhysicalConstants.c
+    h = PhysicalConstants.h
+    ħ = PhysicalConstants.ħ
+    μ0 = PhysicalConstants.μ0
+    ϵ0 = PhysicalConstants.ϵ0
+
+    @test h / ħ ≈ 2 * π
+    @test μ0 / (4.0e-7 * π) ≈ 1.0
+    @test c^2 * μ0 * ϵ0 ≈ 1.0
+
+    @testset "convert unit" begin
+        V = 100 * rand(Float64)
+        _unit_list = [:J, :eV, :meV, :MHz, :GHz, :K, :mK]
+        for origin in _unit_list
+            for middle in _unit_list
+                for target in _unit_list
+                    V_middle = convert_unit(V, origin, middle)
+                    V_target = convert_unit(V_middle, middle, target)
+                    V_origin = convert_unit(V_target, target, origin)
+                    @test V ≈ V_origin
+                end
+            end
+        end
+        @test_throws ArgumentError convert_unit(V, :bad_unit, :J)
+        @test_throws ArgumentError convert_unit(V, :J, :bad_unit)
+    end
+
+    @testset "Type Inference" begin
+        v = rand(Float64)
+        _unit_list = [:J, :eV, :meV, :GHz, :mK]
+        for u1 in _unit_list
+            for u2 in _unit_list
+                @inferred convert_unit(v, u1, u2)
+            end
+        end
+    end
+end
diff --git a/lib/QuantumToolboxUtils/test/cite.jl b/lib/QuantumToolboxUtils/test/cite.jl
new file mode 100644
index 000000000..82a7873b8
--- /dev/null
+++ b/lib/QuantumToolboxUtils/test/cite.jl
@@ -0,0 +1,21 @@
+@testitem "cite" begin
+
+    # citation bibtex
+    io_buffer = IOBuffer()
+    QuantumToolboxUtils.cite(io_buffer)
+    captured_output = String(take!(io_buffer))
+    @test captured_output ==
+        """@article{QuantumToolbox.jl2025,\n""" *
+        """  title = {Quantum{T}oolbox.jl: {A}n efficient {J}ulia framework for simulating open quantum systems},\n""" *
+        """  author = {Mercurio, Alberto and Huang, Yi-Te and Cai, Li-Xun and Chen, Yueh-Nan and Savona, Vincenzo and Nori, Franco},\n""" *
+        """  journal = {{Quantum}},\n""" *
+        """  issn = {2521-327X},\n""" *
+        """  publisher = {{Verein zur F{\\"{o}}rderung des Open Access Publizierens in den Quantenwissenschaften}},\n""" *
+        """  volume = {9},\n""" *
+        """  pages = {1866},\n""" *
+        """  month = sep,\n""" *
+        """  year = {2025},\n""" *
+        """  doi = {10.22331/q-2025-09-29-1866},\n""" *
+        """  url = {https://doi.org/10.22331/q-2025-09-29-1866}\n""" *
+        """}\n"""
+end
diff --git a/lib/QuantumToolboxUtils/test/code_quality.jl b/lib/QuantumToolboxUtils/test/code_quality.jl
new file mode 100644
index 000000000..c31a5cb8f
--- /dev/null
+++ b/lib/QuantumToolboxUtils/test/code_quality.jl
@@ -0,0 +1,9 @@
+@testset "Code quality (QuantumToolboxUtils)" verbose = true begin
+    @testset "Aqua.jl" begin
+        Aqua.test_all(QuantumToolboxUtils; ambiguities = false, unbound_args = false)
+    end
+
+    @testset "JET.jl" begin
+        JET.test_package(QuantumToolboxUtils; target_modules = (QuantumToolboxUtils,), ignore_missing_comparison = true)
+    end
+end
diff --git a/lib/QuantumToolboxUtils/test/linalg.jl b/lib/QuantumToolboxUtils/test/linalg.jl
new file mode 100644
index 000000000..31f412628
--- /dev/null
+++ b/lib/QuantumToolboxUtils/test/linalg.jl
@@ -0,0 +1,38 @@
+@testitem "linalg helpers" begin
+    # meshgrid
+    x = [1, 2, 3]
+    y = [10, 20]
+    X, Y = meshgrid(x, y)
+    @test X == [1 2 3; 1 2 3]
+    @test Y == [10 10 10; 20 20 20]
+end
+
+@testitem "Arnoldi and expv" begin
+    import LinearAlgebra: norm
+
+    # Use n = 3 and m = 2 to avoid exact Arnoldi breakdown (and NaN in AS).
+    A = 0.01 * [1.0 2.0 0.0; 0.0 3.0 4.0; 5.0 0.0 6.0]
+    b = [1.0, 2.0, 3.0]
+    m = 2
+
+    AS = arnoldi(A, b, m)
+    @test AS.m == m
+    @test size(AS.V) == (3, 3)
+    @test size(AS.H) == (3, 2)
+    @test norm(AS.V[:, 1]) ≈ 1.0
+
+    x = expv(A, 0.2, b; m = m)
+    @test x ≈ exp(0.2 * A) * b atol = 1.0e-4
+
+    x2 = zeros(3)
+    expv!(x2, AS, 0.2, b)
+    @test x2 ≈ exp(0.2 * A) * b atol = 1.0e-4
+
+    # test dimension mismatch
+    A_bad = [1.0 0.0; 0.0 1.0]
+    b_bad = [1.0, 0.0]
+    AS_bad_V = ArnoldiSpace(zeros(3, 3), zeros(3, 2), zeros(3, 2), 2)
+    AS_bad_b = ArnoldiSpace(zeros(2, 3), zeros(3, 2), zeros(3, 2), 2)
+    @test_throws DimensionMismatch arnoldi!(AS_bad_V, A_bad, b_bad)
+    @test_throws DimensionMismatch arnoldi!(AS_bad_b, A_bad, [1.0, 0.0, 0.0])
+end
diff --git a/lib/QuantumToolboxUtils/test/physics_func.jl b/lib/QuantumToolboxUtils/test/physics_func.jl
new file mode 100644
index 000000000..6a8fb879e
--- /dev/null
+++ b/lib/QuantumToolboxUtils/test/physics_func.jl
@@ -0,0 +1,33 @@
+@testitem "gaussian" begin
+    mu = 0.7
+    sigma = 1.2
+    delta = 0.4
+
+    @test gaussian(mu, mu, sigma) == 1.0
+    @test gaussian(mu + delta, mu, sigma) ≈ gaussian(mu - delta, mu, sigma)
+
+    @testset "Type Inference" begin
+        @inferred gaussian(0.1, 0.2, 0.3)
+    end
+end
+
+@testitem "n_thermal" begin
+
+    ω1 = rand(Float64)
+    ω2 = rand(Float64)
+    @test n_thermal(0, ω2) == 0.0
+    @test n_thermal(ω1, 0) == 0.0
+    @test n_thermal(ω1, -ω2) == 0.0
+    @test n_thermal(ω1, ω2) == 1 / (exp(ω1 / ω2) - 1)
+    @test typeof(n_thermal(Int32(2), Int32(3))) == Float32
+    @test typeof(n_thermal(Float32(2), Float32(3))) == Float32
+    @test typeof(n_thermal(Int64(2), Int32(3))) == Float64
+    @test typeof(n_thermal(Int32(2), Int64(3))) == Float64
+    @test typeof(n_thermal(Float64(2), Float32(3))) == Float64
+    @test typeof(n_thermal(Float32(2), Float64(3))) == Float64
+
+    @testset "Type Inference" begin
+        v = rand(Float64)
+        @inferred n_thermal(v, Int32(123))
+    end
+end
diff --git a/lib/QuantumToolboxUtils/test/versioninfo_and_settings.jl b/lib/QuantumToolboxUtils/test/versioninfo_and_settings.jl
new file mode 100644
index 000000000..36578cbcb
--- /dev/null
+++ b/lib/QuantumToolboxUtils/test/versioninfo_and_settings.jl
@@ -0,0 +1,33 @@
+@testitem "versioninfo and about" begin
+    # test error handling
+    io_wrong = IOBuffer()
+    push!(QuantumToolboxUtils._QT_LIBRARIES, Base)
+    @test_throws ArgumentError QuantumToolboxUtils.versioninfo(io_wrong)
+    deleteat!(QuantumToolboxUtils._QT_LIBRARIES, findall(x -> x == Base, QuantumToolboxUtils._QT_LIBRARIES))
+
+    # versioninfo
+    io_version = IOBuffer()
+    QuantumToolboxUtils.versioninfo(io_version)
+    version_output = String(take!(io_version))
+    @test occursin("QuantumToolbox.jl: Quantum Toolbox in Julia", version_output)
+    @test occursin("Package information:", version_output)
+    @test occursin("QuantumToolboxUtils", version_output)
+    @test occursin("System information:", version_output)
+    @test occursin("Please cite QuantumToolbox.jl in your publication", version_output)
+
+    # about
+    io_about = IOBuffer()
+    QuantumToolboxUtils.about(io_about)
+    about_output = String(take!(io_about))
+    @test version_output == about_output
+end
+
+@testitem "Settings" begin
+    io = IOBuffer()
+    show(io, QuantumToolboxUtils.settings)
+    out = String(take!(io))
+    @test occursin("QuantumToolbox.jl Settings", out)
+    @test occursin("tidyup_tol", out)
+    @test occursin("auto_tidyup", out)
+    @test occursin("ProgressMeterKWARGS", out)
+end
diff --git a/src/QuantumToolbox.jl b/src/QuantumToolbox.jl
index e4a0a740c..1019e2676 100644
--- a/src/QuantumToolbox.jl
+++ b/src/QuantumToolbox.jl
@@ -11,8 +11,29 @@ import Pkg
 import Random: AbstractRNG, default_rng, seed!
 import Statistics: mean, std
 
+## Re-export of QuantumToolbox libraries
+import Reexport: @reexport
+@reexport using QuantumToolboxUtils
+
+## internal functions of QuantumToolbox libraries
+import QuantumToolboxUtils:
+    FloatOrComplex,
+    getVal,
+    makeVal,
+    get_typename_wrapper,
+    _float_type,
+    _complex_float_type,
+    _convert_eltype_wordsize,
+    _non_static_array_warning,
+    _lazy_tensor_warning,
+    _Ginibre_ensemble,
+    _Boltzmann_weight,
+    _dense_similar,
+    _sparse_similar
+
 ## SciML packages (for QobjEvo, OrdinaryDiffEq, and LinearSolve)
 import SciMLBase:
+    SciMLBase,
     solve,
     solve!,
     init,
@@ -87,11 +108,8 @@ export cache_operator, iscached, isconstant
 
 # Source files
 
-## Utility
-include("settings.jl")
-include("utilities.jl")
-include("versioninfo.jl")
-include("linear_maps.jl")
+## Some overloading with QuantumToolboxUtils library
+include("utils.jl")
 
 ## Quantum Object
 include("qobj/dimensions.jl")
@@ -131,7 +149,6 @@ include("time_evolution/time_evolution_dynamical.jl")
 include("correlations.jl")
 include("wigner.jl")
 include("spin_lattice.jl")
-include("arnoldi.jl")
 include("entropy.jl")
 include("metrics.jl")
 include("negativity.jl")
@@ -147,4 +164,9 @@ include("visualization/wigner.jl")
 ## deprecated functions
 include("deprecated.jl")
 
+function __init__()
+    QuantumToolboxUtils._register_qt_library!(QuantumToolbox)
+    return nothing
+end
+
 end
diff --git a/src/utilities.jl b/src/utilities.jl
deleted file mode 100644
index 0a58f429b..000000000
--- a/src/utilities.jl
+++ /dev/null
@@ -1,204 +0,0 @@
-#=
-Utilities:
-    internal (or external) functions which will be used throughout the entire package
-=#
-
-export gaussian, n_thermal
-export PhysicalConstants, convert_unit
-export row_major_reshape, meshgrid
-
-@doc raw"""
-    row_major_reshape(Q::AbstractArray, shapes...)
-
-Reshapes `Q` in the row-major order, as numpy.
-"""
-row_major_reshape(Q::AbstractArray{T}, shapes...) where {T} =
-    PermutedDimsArray(reshape(Q, reverse(shapes)...), (length(shapes):-1:1))
-
-@doc raw"""
-    meshgrid(x::AbstractVector, y::AbstractVector)
-
-Equivalent to [numpy meshgrid](https://numpy.org/doc/stable/reference/generated/numpy.meshgrid.html).
-"""
-function meshgrid(x::AbstractVector{T}, y::AbstractVector{T}) where {T}
-    X = reshape(repeat(x, inner = length(y)), length(y), length(x))
-    Y = repeat(y, outer = (1, length(x)))
-    return X, Y
-end
-
-@doc raw"""
-    gaussian(x::Number, μ::Number, σ::Number)
-
-Returns the gaussian function ``\exp \left[- \frac{(x - \mu)^2}{2 \sigma^2} \right]``,
-where ``\mu`` and ``\sigma^2`` are the mean and the variance respectively.
-"""
-gaussian(x::Number, μ::Number, σ::Number) = exp(-(x - μ)^2 / (2 * σ^2))
-
-@doc raw"""
-    n_thermal(ω::Real, ω_th::Real)
-
-Return the number of photons in thermal equilibrium for an harmonic oscillator mode with frequency ``\omega``, at the temperature described by ``\omega_{\textrm{th}} \equiv k_B T / \hbar``:
-```math
-n(\omega, \omega_{\textrm{th}}) = \frac{1}{e^{\omega/\omega_{\textrm{th}}} - 1},
-```
-where ``\hbar`` is the reduced Planck constant, and ``k_B`` is the Boltzmann constant.
-"""
-function n_thermal(ω::T1, ω_th::T2) where {T1 <: Real, T2 <: Real}
-    x = exp(ω / ω_th)
-    n = ((x != 1) && (ω_th > 0)) ? 1 / (x - 1) : 0
-    return _float_type(promote_type(T1, T2))(n)
-end
-
-@doc raw"""
-    const PhysicalConstants
-
-A `NamedTuple` which stores some constant values listed in [*CODATA recommended values of the fundamental physical constants: 2022*](https://physics.nist.gov/cuu/pdf/wall_2022.pdf)
-
-The current stored constants are:
-- `c` : (exact) speed of light in vacuum with unit ``[\textrm{m}\cdot\textrm{s}^{-1}]``
-- `G` : Newtonian constant of gravitation with unit ``[\textrm{m}^3\cdot\textrm{kg}^{−1}\cdot\textrm{s}^{−2}]``
-- `h` : (exact) Planck constant with unit ``[\textrm{J}\cdot\textrm{s}]``
-- `ħ` : reduced Planck constant with unit ``[\textrm{J}\cdot\textrm{s}]``
-- `e` : (exact) elementary charge with unit ``[\textrm{C}]``
-- `μ0` : vacuum magnetic permeability with unit ``[\textrm{N}\cdot\textrm{A}^{-2}]``
-- `ϵ0` : vacuum electric permittivity with unit ``[\textrm{F}\cdot\textrm{m}^{-1}]``
-- `k` : (exact) Boltzmann constant with unit ``[\textrm{J}\cdot\textrm{K}^{-1}]``
-- `NA` : (exact) Avogadro constant with unit ``[\textrm{mol}^{-1}]``
-
-# Examples
-
-```jldoctest
-julia> PhysicalConstants.ħ
-1.0545718176461565e-34
-```
-"""
-const PhysicalConstants = (
-    c = 299792458.0,
-    G = 6.6743e-11,
-    h = 6.62607015e-34,
-    ħ = 6.62607015e-34 / (2 * π),
-    e = 1.602176634e-19,
-    μ0 = 1.25663706127e-6,
-    ϵ0 = 8.8541878188e-12,
-    k = 1.380649e-23,
-    NA = 6.02214076e23,
-)
-
-# common energy units (the values below are all in the unit of Joule)
-const _energy_units::Dict{Symbol, Float64} = Dict(
-    :J => 1.0,
-    :eV => PhysicalConstants.e,
-    :meV => 1.0e-3 * PhysicalConstants.e,
-    :MHz => 1.0e6 * PhysicalConstants.h,
-    :GHz => 1.0e9 * PhysicalConstants.h,
-    :K => PhysicalConstants.k,
-    :mK => 1.0e-3 * PhysicalConstants.k,
-)
-
-@doc raw"""
-    convert_unit(value::Real, unit1::Symbol, unit2::Symbol)
-
-Convert the energy `value` from `unit1` to `unit2`. The `unit1` and `unit2` can be either the following `Symbol`:
-- `:J` : Joule
-- `:eV` : electron volt
-- `:meV` : milli-electron volt
-- `:MHz` : Mega-Hertz multiplied by Planck constant ``h``
-- `:GHz` : Giga-Hertz multiplied by Planck constant ``h``
-- `:K` : Kelvin multiplied by Boltzmann constant ``k``
-- `:mK` : milli-Kelvin multiplied by Boltzmann constant ``k``
-
-Note that we use the values stored in [`PhysicalConstants`](@ref) to do the conversion.
-
-# Examples
-
-```jldoctest
-julia> convert_unit(1, :eV, :J)
-1.602176634e-19
-
-julia> convert_unit(1, :GHz, :J)
-6.62607015e-25
-
-julia> round(convert_unit(1, :meV, :mK), digits=4)
-11604.5181
-```
-"""
-function convert_unit(value::T, unit1::Symbol, unit2::Symbol) where {T <: Real}
-    !haskey(_energy_units, unit1) && throw(ArgumentError("Invalid unit :$(unit1)"))
-    !haskey(_energy_units, unit2) && throw(ArgumentError("Invalid unit :$(unit2)"))
-    return _float_type(T)(value * (_energy_units[unit1] / _energy_units[unit2]))
-end
-
-get_typename_wrapper(A) = Base.typename(typeof(A)).wrapper
-
-_dense_similar(A::AbstractArray, args...) = similar(A, args...)
-_dense_similar(A::AbstractSparseMatrix, args...) = similar(nonzeros(A), args...)
-
-_sparse_similar(A::AbstractArray, args...) = sparse(args...)
-
-_Ginibre_ensemble(rng::AbstractRNG, ::Type{T}, n::Int, rank::Int = n) where {T <: Complex} = randn(rng, T, n, rank) / sqrt(T(n))
-
-_Boltzmann_weight(β::T, E::Int) where {T <: Real} = (E != 0 || isfinite(β)) ? exp(-β * E) : one(T)
-
-makeVal(x::Val{T}) where {T} = x
-makeVal(x) = Val(x)
-
-getVal(x::Val{T}) where {T} = T
-getVal(x) = x # getVal for any other type
-
-_non_static_array_warning(argname, arg::Tuple{}) =
-    throw(ArgumentError("The argument $argname must be a Tuple or a StaticVector of non-zero length."))
-_non_static_array_warning(argname, arg::Union{SVector{N, T}, MVector{N, T}, NTuple{N, T}}) where {N, T} = nothing
-_non_static_array_warning(argname, arg::AbstractVector{T}) where {T} =
-    @warn "The argument $argname should be a Tuple or a StaticVector for better performance. Try to use `$argname = $(Tuple(arg))` instead of `$argname = $arg`. " *
-    "Alternatively, you can do `import QuantumToolbox: SVector` " *
-    "and use `$argname = SVector(" *
-    join(arg, ", ") *
-    ")`." maxlog = 1
-
-# lazy tensor warning
-for AType in (:AbstractArray, :AbstractSciMLOperator)
-    for BType in (:AbstractArray, :AbstractSciMLOperator)
-        if AType == BType == :AbstractArray
-            @eval begin
-                _lazy_tensor_warning(::$AType, ::$BType) = nothing
-            end
-        else
-            @eval begin
-                _lazy_tensor_warning(A::$AType, B::$BType) =
-                    @warn "using lazy tensor (which can hurt performance) between data types: $(get_typename_wrapper(A)) and $(get_typename_wrapper(B))"
-            end
-        end
-    end
-end
-
-# alias of abstract types
-const FloatOrComplex = Union{AbstractFloat, Complex}
-
-# functions for getting Float or Complex element type
-_float_type(::AbstractArray{T}) where {T <: Number} = _float_type(T)
-_float_type(::AbstractSciMLOperator{T}) where {T <: Number} = _float_type(T)
-_float_type(::Type{Int32}) = Float32
-_float_type(::Type{Int64}) = Float64
-_float_type(::Type{Complex{Int32}}) = Float32
-_float_type(::Type{Complex{Int64}}) = Float64
-_float_type(::Type{Complex{T}}) where {T <: Real} = T
-_float_type(T::Type{<:AbstractFloat}) = T # Allow other untracked Real types, like ForwardDiff.Dual
-_complex_float_type(::AbstractArray{T}) where {T <: Number} = _complex_float_type(T)
-_complex_float_type(::AbstractSciMLOperator{T}) where {T <: Number} = _complex_float_type(T)
-_complex_float_type(::Type{Int32}) = ComplexF32
-_complex_float_type(::Type{Int64}) = ComplexF64
-_complex_float_type(::Type{Float32}) = ComplexF32
-_complex_float_type(::Type{Float64}) = ComplexF64
-_complex_float_type(::Type{Complex{Int32}}) = ComplexF32
-_complex_float_type(::Type{Complex{Int64}}) = ComplexF64
-_complex_float_type(::Type{Complex{Float32}}) = ComplexF32
-_complex_float_type(::Type{Complex{Float64}}) = ComplexF64
-_complex_float_type(T::Type{<:AbstractFloat}) = Complex{T} # Allow other untracked Complex types, like ForwardDiff.Dual
-_complex_float_type(T::Type{<:Complex}) = T       # Allow other untracked Complex types, like ForwardDiff.Dual
-
-_convert_eltype_wordsize(::Type{T}, ::Val{64}) where {T <: Int} = Int64
-_convert_eltype_wordsize(::Type{T}, ::Val{32}) where {T <: Int} = Int32
-_convert_eltype_wordsize(::Type{T}, ::Val{64}) where {T <: AbstractFloat} = Float64
-_convert_eltype_wordsize(::Type{T}, ::Val{32}) where {T <: AbstractFloat} = Float32
-_convert_eltype_wordsize(::Type{Complex{T}}, ::Val{64}) where {T <: Union{Int, AbstractFloat}} = ComplexF64
-_convert_eltype_wordsize(::Type{Complex{T}}, ::Val{32}) where {T <: Union{Int, AbstractFloat}} = ComplexF32
diff --git a/src/utils.jl b/src/utils.jl
new file mode 100644
index 000000000..a9236d4fc
--- /dev/null
+++ b/src/utils.jl
@@ -0,0 +1,13 @@
+#=
+Some alias and overloading functions from QuantumToolboxUtils
+=#
+
+function QuantumToolboxUtils._add_library_deps!(lib::Val{:QuantumToolbox}, DEPpkgs::Vector{Module})
+    QuantumToolboxUtils._add_pkgs!(DEPpkgs, Module[SciMLBase, SciMLOperators, OrdinaryDiffEqCore, LinearSolve])
+    return nothing
+end
+
+const settings = QuantumToolboxUtils.settings
+const versioninfo = QuantumToolboxUtils.versioninfo
+const about = QuantumToolboxUtils.about
+const cite = QuantumToolboxUtils.cite
diff --git a/test/core-test/code-quality/Project.toml b/test/code-quality/Project.toml
similarity index 64%
rename from test/core-test/code-quality/Project.toml
rename to test/code-quality/Project.toml
index d71385453..5a4113b71 100644
--- a/test/core-test/code-quality/Project.toml
+++ b/test/code-quality/Project.toml
@@ -2,8 +2,12 @@
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
 QuantumToolbox = "6c2fb7c5-b903-41d2-bc5e-5a7c320b9fab"
+QuantumToolboxUtils = "02afad9e-e2b8-4f9a-965f-b5dca5b538ca"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
+[sources]
+QuantumToolboxUtils = {path = "../../lib/QuantumToolboxUtils"}
+
 [compat]
 Aqua = "0.8"
 JET = "0.9, 0.10, 0.11"
diff --git a/test/core-test/code-quality/code_quality.jl b/test/code-quality/code_quality.jl
similarity index 81%
rename from test/core-test/code-quality/code_quality.jl
rename to test/code-quality/code_quality.jl
index 71df061d4..f2d416f14 100644
--- a/test/core-test/code-quality/code_quality.jl
+++ b/test/code-quality/code_quality.jl
@@ -1,4 +1,4 @@
-@testset "Code quality" verbose = true begin
+@testset "Code quality (QuantumToolbox)" verbose = true begin
     @testset "Aqua.jl" begin
         Aqua.test_all(QuantumToolbox; ambiguities = false, unbound_args = false)
     end
diff --git a/test/core-test/utilities.jl b/test/core-test/utilities.jl
deleted file mode 100644
index e0cdf96e1..000000000
--- a/test/core-test/utilities.jl
+++ /dev/null
@@ -1,77 +0,0 @@
-@testitem "Utilities" begin
-
-    # citation bibtex
-    io_buffer = IOBuffer()
-    QuantumToolbox.cite(io_buffer)
-    captured_output = String(take!(io_buffer))
-    @test captured_output ==
-        """@article{QuantumToolbox.jl2025,\n""" *
-        """  title = {Quantum{T}oolbox.jl: {A}n efficient {J}ulia framework for simulating open quantum systems},\n""" *
-        """  author = {Mercurio, Alberto and Huang, Yi-Te and Cai, Li-Xun and Chen, Yueh-Nan and Savona, Vincenzo and Nori, Franco},\n""" *
-        """  journal = {{Quantum}},\n""" *
-        """  issn = {2521-327X},\n""" *
-        """  publisher = {{Verein zur F{\\"{o}}rderung des Open Access Publizierens in den Quantenwissenschaften}},\n""" *
-        """  volume = {9},\n""" *
-        """  pages = {1866},\n""" *
-        """  month = sep,\n""" *
-        """  year = {2025},\n""" *
-        """  doi = {10.22331/q-2025-09-29-1866},\n""" *
-        """  url = {https://doi.org/10.22331/q-2025-09-29-1866}\n""" *
-        """}\n"""
-
-    @testset "n_thermal" begin
-        ω1 = rand(Float64)
-        ω2 = rand(Float64)
-        @test n_thermal(0, ω2) == 0.0
-        @test n_thermal(ω1, 0) == 0.0
-        @test n_thermal(ω1, -ω2) == 0.0
-        @test n_thermal(ω1, ω2) == 1 / (exp(ω1 / ω2) - 1)
-        @test typeof(n_thermal(Int32(2), Int32(3))) == Float32
-        @test typeof(n_thermal(Float32(2), Float32(3))) == Float32
-        @test typeof(n_thermal(Int64(2), Int32(3))) == Float64
-        @test typeof(n_thermal(Int32(2), Int64(3))) == Float64
-        @test typeof(n_thermal(Float64(2), Float32(3))) == Float64
-        @test typeof(n_thermal(Float32(2), Float64(3))) == Float64
-    end
-
-    @testset "CODATA Physical Constants" begin
-        c = PhysicalConstants.c
-        h = PhysicalConstants.h
-        ħ = PhysicalConstants.ħ
-        μ0 = PhysicalConstants.μ0
-        ϵ0 = PhysicalConstants.ϵ0
-
-        @test h / ħ ≈ 2 * π
-        @test μ0 / (4.0e-7 * π) ≈ 1.0
-        @test c^2 * μ0 * ϵ0 ≈ 1.0
-    end
-
-    @testset "convert unit" begin
-        V = 100 * rand(Float64)
-        _unit_list = [:J, :eV, :meV, :MHz, :GHz, :K, :mK]
-        for origin in _unit_list
-            for middle in _unit_list
-                for target in _unit_list
-                    V_middle = convert_unit(V, origin, middle)
-                    V_target = convert_unit(V_middle, middle, target)
-                    V_origin = convert_unit(V_target, target, origin)
-                    @test V ≈ V_origin
-                end
-            end
-        end
-        @test_throws ArgumentError convert_unit(V, :bad_unit, :J)
-        @test_throws ArgumentError convert_unit(V, :J, :bad_unit)
-    end
-
-    @testset "Type Inference" begin
-        v1 = rand(Float64)
-        @inferred n_thermal(v1, Int32(123))
-
-        _unit_list = [:J, :eV, :meV, :GHz, :mK]
-        for u1 in _unit_list
-            for u2 in _unit_list
-                @inferred convert_unit(v1, u1, u2)
-            end
-        end
-    end
-end
diff --git a/test/runtests.jl b/test/runtests.jl
index 3cad44db2..dd3795ea6 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -4,21 +4,28 @@ using Pkg
 
 const testdir = dirname(@__FILE__)
 
+# Define the paths to the library
+const LIBRARY_NAME_AND_PATH = Dict(
+    "Utils" => ("QuantumToolboxUtils", joinpath(testdir, "..", "lib", "QuantumToolboxUtils")),
+)
+const LIBRARY_LIST = collect(keys(LIBRARY_NAME_AND_PATH))
+
 # Define the paths to the extension tests
-const EXTENSION_PATHS = Dict(
-    "AutoDiff_Ext" => joinpath(testdir, "ext-test", "cpu", "autodiff"),
-    "Makie_Ext" => joinpath(testdir, "ext-test", "cpu", "makie"),
-    "CUDA_Ext" => joinpath(testdir, "ext-test", "gpu"),
+const EXTENSION_PATH = Dict(
+    "AutoDiff-Ext" => joinpath(testdir, "ext-test", "cpu", "autodiff"),
+    "Makie-Ext" => joinpath(testdir, "ext-test", "cpu", "makie"),
+    "CUDA-Ext" => joinpath(testdir, "ext-test", "gpu"),
     "Arbitrary-Precision" => joinpath(testdir, "ext-test", "cpu", "arbitrary_precision"),
 )
-const EXTENSION_LIST = collect(keys(EXTENSION_PATHS))
+const EXTENSION_LIST = collect(keys(EXTENSION_PATH))
 
 # Handle the GROUP environment variable to determine which tests to run
 const GROUP = get(ENV, "GROUP", "All")
 const GROUP_LIST = String[
     "All",
-    "Core",
+    "Main",
     "Code-Quality",
+    LIBRARY_LIST...,
     EXTENSION_LIST...,
 ]
 (GROUP in GROUP_LIST) || throw(ArgumentError("Unknown GROUP = $GROUP\nThe allowed groups are: $GROUP_LIST\n"))
@@ -26,21 +33,51 @@ const GROUP_LIST = String[
 # function to set up the environment for subtests
 function setup_subtest_env(path::String)
     Pkg.activate(path)
+    if VERSION < v"1.11"
+        for lib in LIBRARY_LIST
+            _, lib_path = LIBRARY_NAME_AND_PATH[lib]
+            Pkg.develop(path = lib_path)
+        end
+    end
     Pkg.develop(PackageSpec(path = dirname(@__DIR__)))
     Pkg.update()
     return nothing
 end
 
-######################################
-# Core tests (use TestItemRunner.jl) #
-######################################
-if (GROUP == "All") || (GROUP == "Core")
+#######################################
+# Main tests (use TestItemRunner.jl) #
+#######################################
+if (GROUP == "All") || (GROUP == "Main")
+    import QuantumToolbox
+
+    QuantumToolbox.about()
+
+    println("\nStart running Main tests...")
+
+    # tests in lib folder for each library
+    # PATH: lib/LIBRARY_NAME/test/
+    for lib in LIBRARY_LIST
+        lib_name, _ = LIBRARY_NAME_AND_PATH[lib]
+        println("\n[$lib_name]")
+        @run_package_tests filter = ti -> occursin(joinpath("lib", lib_name, "test"), ti.filename) verbose = true
+    end
+
+    # main package tests (all tests except those in the lib folder)
+    println("\n[QuantumToolbox]")
+    @run_package_tests filter = ti -> !occursin("lib", ti.filename) verbose = true
+
+    println("\n===============> Main tests completed <===============\n")
+end
+
+# only run tests for a specific library
+if GROUP ∈ LIBRARY_LIST
     import QuantumToolbox
 
     QuantumToolbox.about()
 
-    println("\nStart running Core tests...\n")
-    @run_package_tests verbose = true
+    lib_name, _ = LIBRARY_NAME_AND_PATH[GROUP]
+    println("\n[$lib_name]")
+    @run_package_tests filter = ti -> occursin(joinpath("lib", lib_name, "test"), ti.filename) verbose = true
 end
 
 ########################################################################
@@ -49,7 +86,7 @@ end
 
 # Code Quality tests
 if (GROUP == "All") || (GROUP == "Code-Quality")
-    path = joinpath(testdir, "core-test", "code-quality")
+    path = joinpath(testdir, "code-quality")
     setup_subtest_env(path)
 
     using QuantumToolbox
@@ -57,20 +94,25 @@ if (GROUP == "All") || (GROUP == "Code-Quality")
 
     (GROUP == "Code-Quality") && QuantumToolbox.about() # print version info. for code quality CI in GitHub
 
+    # run code quality tests for all libraries and the main package
+    for lib in LIBRARY_LIST
+        _, lib_path = LIBRARY_NAME_AND_PATH[lib]
+        include(joinpath(lib_path, "test", "code_quality.jl"))
+    end
     include(joinpath(path, "code_quality.jl"))
 end
 
 # Extension tests
 if GROUP ∈ EXTENSION_LIST
-    path = EXTENSION_PATHS[GROUP]
+    path = EXTENSION_PATH[GROUP]
     setup_subtest_env(path)
 
-    if GROUP == "AutoDiff_Ext"
+    if GROUP == "AutoDiff-Ext"
         println(Pkg.status())
         include(joinpath(path, "autodiff.jl"))
-    elseif GROUP == "Makie_Ext"
+    elseif GROUP == "Makie-Ext"
         include(joinpath(path, "makie_ext.jl"))
-    elseif GROUP == "CUDA_Ext"
+    elseif GROUP == "CUDA-Ext"
         include(joinpath(path, "cuda_ext.jl"))
     elseif GROUP == "Arbitrary-Precision"
         include(joinpath(path, "arbitrary_precision.jl"))