diff --git a/src/QuantumStateTransfer.jl b/src/QuantumStateTransfer.jl
index 7e396e3..5e2c459 100644
--- a/src/QuantumStateTransfer.jl
+++ b/src/QuantumStateTransfer.jl
@@ -22,9 +22,9 @@ include("types.jl")
 include("EpsilonOptimization/EpsilonOptimization.jl")
 using .EpsilonOptimization
 
-include("state_transfer.jl")
-include("uniform_mixing.jl")
-include("fractional_revival.jl")
+include("core/state_transfer.jl")
+include("core/uniform_mixing.jl")
+include("core/fractional_revival.jl")
 
 # TODO: Exports (add more later)
 export max_state_transfer, check_state_transfer
diff --git a/src/fractional_revival.jl b/src/core/fractional_revival.jl
similarity index 100%
rename from src/fractional_revival.jl
rename to src/core/fractional_revival.jl
diff --git a/src/state_transfer.jl b/src/core/state_transfer.jl
similarity index 96%
rename from src/state_transfer.jl
rename to src/core/state_transfer.jl
index b616831..c26c107 100644
--- a/src/state_transfer.jl
+++ b/src/core/state_transfer.jl
@@ -131,8 +131,7 @@ end
 [TODO: Write here]
 
 # Notes
-[TODO: Write here. Proof sketch of bounds on Lipschitz and alpha constants, plus further
-relevant references?]
+[TODO: Refer to [`transfer_fidelity_deriv_bound`](@ref) for proof sketch of bounds]
 """
 function max_state_transfer(g::AbstractGraph, args...)
     if !is_simple(g)
@@ -151,7 +150,7 @@ function max_state_transfer(
     epsilon::Real,
     method::Symbol=:lipschitz_bb,
 ) where {Tl<:Union{Integer,Tuple{Integer,Integer}}}
-    if !is_zero_diagonal_symmetric(A)
+    if !is_zero_diag_symmetric(A)
         throw(ArgumentError("Matrix must be symmetric with zero diagonal"))
     end
 
@@ -194,8 +193,7 @@ end
 [TODO: Write here]
 
 # Notes
-[TODO: Write here. Proof sketch of bounds on Lipschitz and alpha constants, plus further
-relevant references?]
+[TODO: Refer to [`transfer_fidelity_deriv_bound`](@ref) for proof sketch of bounds]
 """
 function check_state_transfer(g::AbstractGraph, args...)
     if !is_simple(g)
@@ -215,7 +213,7 @@ function check_state_transfer(
     epsilon::Real,
     method::Symbol=:lipschitz_bb,
 ) where {Tl<:Union{Integer,Tuple{Integer,Integer}}}
-    if !is_zero_diagonal_symmetric(A)
+    if !is_zero_diag_symmetric(A)
         throw(ArgumentError("Matrix must be symmetric with zero diagonal"))
     end
 
@@ -352,12 +350,12 @@ function _optimize_state_transfer_impl(input::_StateTransferProblemInput)
     end
 
     if method == :lipschitz_bb
-        lipschitz_constant = maximum(norm.(eachcol(A)))
+        lipschitz_constant = transfer_fidelity_deriv_bound(A, 1)
         solver = LipschitzBranchAndBound(
             epsilon, lipschitz_constant; target=target_infidelity
         )
     elseif method == :alpha_bb
-        alpha = maximum(norm.(eachcol(A^2))) / 2
+        alpha = transfer_fidelity_deriv_bound(A, 2) / 2
         solver = AlphaBranchAndBound(epsilon, alpha; target=target_infidelity)
     else
         throw(
diff --git a/src/uniform_mixing.jl b/src/core/uniform_mixing.jl
similarity index 100%
rename from src/uniform_mixing.jl
rename to src/core/uniform_mixing.jl
diff --git a/src/utils.jl b/src/utils.jl
index 4ea82a3..37824fd 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -5,7 +5,7 @@
 # distributed except according to those terms.
 
 """
-    is_zero_diagonal_symmetric(A) -> Bool
+    is_zero_diag_symmetric(A) -> Bool
 
 Check whether a matrix `A` is symmetric with a zero diagonal.
 
@@ -24,7 +24,7 @@ only consider here real-valued adjacency matrices) is required for the walk Hami
 # Examples
 [TODO: Write here]
 """
-function is_zero_diagonal_symmetric(A::AbstractMatrix{<:Real})
+function is_zero_diag_symmetric(A::AbstractMatrix{<:Real})
     (m, n) = size(A)
 
     return m == n && # Square
@@ -54,3 +54,21 @@ adjacency matrix of `g`.
 function is_simple(g::AbstractGraph)
     return !isdirected(g) && !has_self_loops(g)
 end
+
+"""
+    transfer_fidelity_deriv_bound(A, order) -> Float64
+
+[TODO: Write here]
+
+# Arguments
+[TODO: Write here]
+
+# Returns
+[TODO: Write here]
+
+# Notes
+[TODO: Proof sketch of bound, plus further relevant references?]
+"""
+function transfer_fidelity_deriv_bound(A::Matrix{Float64}, order::Int)
+    return maximum(norm.eachcol(A^order))
+end