diff --git a/CHANGELOG.md b/CHANGELOG.md
index 3d4d56ddc..e06ec247b 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -7,7 +7,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased](https://github.com/qutip/QuantumToolbox.jl/tree/main)
 
-- Make CUDA conversion more general using Adapt.jl. ([#437])
+- Make CUDA conversion more general using Adapt.jl. ([#436], [#437])
+- Make the generation of `fock` states non-mutating to support Zygote.jl. ([#438])
 
 ## [v0.29.1]
 Release date: 2025-03-07
@@ -189,4 +190,6 @@ Release date: 2024-11-13
 [#423]: https://github.com/qutip/QuantumToolbox.jl/issues/423
 [#428]: https://github.com/qutip/QuantumToolbox.jl/issues/428
 [#430]: https://github.com/qutip/QuantumToolbox.jl/issues/430
+[#436]: https://github.com/qutip/QuantumToolbox.jl/issues/436
 [#437]: https://github.com/qutip/QuantumToolbox.jl/issues/437
+[#438]: https://github.com/qutip/QuantumToolbox.jl/issues/438
diff --git a/src/qobj/states.jl b/src/qobj/states.jl
index 886982f7c..784c172aa 100644
--- a/src/qobj/states.jl
+++ b/src/qobj/states.jl
@@ -37,8 +37,7 @@ function fock(N::Int, j::Int = 0; dims::Union{Int,AbstractVector{Int},Tuple} = N
     if getVal(sparse)
         array = sparsevec([j + 1], [1.0 + 0im], N)
     else
-        array = zeros(ComplexF64, N)
-        array[j+1] = 1
+        array = [i == (j + 1) ? 1.0 + 0im : 0.0 + 0im for i in 1:N]
     end
     return QuantumObject(array; type = Ket, dims = dims)
 end