diff --git a/src/tensor_operations/matrix_algebra.jl b/src/tensor_operations/matrix_algebra.jl
index 4b125503b5..4099c783c5 100644
--- a/src/tensor_operations/matrix_algebra.jl
+++ b/src/tensor_operations/matrix_algebra.jl
@@ -43,46 +43,70 @@ function exp(A::ITensor, Linds, Rinds; ishermitian=false)
     end
   end
 
-  N = ndims(A)
   NL = length(Linds)
   NR = length(Rinds)
   NL != NR && error("Must have equal number of left and right indices")
-  N != NL + NR &&
+  ndims(A) != NL + NR &&
     error("Number of left and right indices must add up to total number of indices")
 
-  # Linds, Rinds may not have the correct directions
-  # TODO: does the need a conversion?
-  Lis = Linds
-  Ris = Rinds
+  # Replace Linds and Rinds with index sets having
+  # same id's but arrow directions as on A
+  Linds = permute(commoninds(A, Linds), Linds)
+  Rinds = permute(commoninds(A, Rinds), Rinds)
 
-  # Ensure the indices have the correct directions,
-  # QNs, etc.
-  # First grab the indices in A, then permute them
-  # correctly.
-  Lis = permute(commoninds(A, Lis), Lis)
-  Ris = permute(commoninds(A, Ris), Ris)
-
-  for (l, r) in zip(Lis, Ris)
+  # Ensure indices have correct directions, QNs, etc.
+  for (l, r) in zip(Linds, Rinds)
     if space(l) != space(r)
       error("In exp, indices must come in pairs with equal spaces.")
     end
-    if hasqns(A)
-      if dir(l) == dir(r)
-        error("In exp, indices must come in pairs with opposite directions")
-      end
+    if hasqns(A) && dir(l) == dir(r)
+      error("In exp, indices must come in pairs with opposite directions")
+    end
+  end
+
+  # <fermions>
+  auto_fermion_enabled = using_auto_fermion()
+  if auto_fermion_enabled
+    # If fermionic, bring indices into i',j',..,dag(j),dag(i)
+    # ordering with Out indices coming before In indices
+    # Resulting tensor acts like a normal matrix (no extra signs
+    # when taking powers A^n)
+    if all(j->dir(j)==Out, Linds) && all(j->dir(j)==In, Rinds)
+      ordered_inds = [Linds..., reverse(Rinds)...]
+    elseif all(j->dir(j)==Out, Rinds) && all(j->dir(j)==In, Linds)
+      ordered_inds = [Rinds..., reverse(Linds)...]
+    else
+      error(
+        "For fermionic exp, Linds and Rinds must have same directions within each set. Got dir.(Linds)=",
+        dir.(Linds),
+        ", dir.(Rinds)=",
+        dir.(Rinds),
+      )
     end
+    A = permute(A, ordered_inds)
+    # A^n now sign free, ok to temporarily disable fermion system
+    disable_auto_fermion()
   end
 
-  CL = combiner(Lis...; dir=Out)
-  CR = combiner(Ris...; dir=In)
+  CL = combiner(Linds...; dir=Out)
+  CR = combiner(Rinds...; dir=In)
   AC = (A * CR) * CL
   expAT = ishermitian ? exp(Hermitian(tensor(AC))) : exp(tensor(AC))
-  return (itensor(expAT) * dag(CR)) * dag(CL)
+  expA = (itensor(expAT) * dag(CR)) * dag(CL)
+
+  # <fermions>
+  if auto_fermion_enabled
+    # Ensure expA indices in "matrix" form before re-enabling fermion system
+    expA = permute(expA, ordered_inds)
+    enable_auto_fermion()
+  end
+
+  return expA
 end
 
 function exp(A::ITensor; kwargs...)
   Ris = filterinds(A; plev=0)
-  Lis = Ris'
+  Lis = dag.(prime.(Ris))
   return exp(A, Lis, Ris; kwargs...)
 end
 
diff --git a/test/base/test_fermions.jl b/test/base/test_fermions.jl
index c70603a283..b68d273663 100644
--- a/test/base/test_fermions.jl
+++ b/test/base/test_fermions.jl
@@ -1,6 +1,6 @@
 using ITensors, Test
 import ITensors: Out, In
-using ITensors.SiteTypes: op, siteinds
+using ITensors.SiteTypes: op, siteind, siteinds
 
 @testset "Fermions" begin
   ITensors.enable_auto_fermion()
@@ -778,5 +778,30 @@ using ITensors.SiteTypes: op, siteinds
     @test norm(a * u - u' * d) ≈ 0 atol = √(eps(real(eltype(a))))
   end
 
+  @testset "Fermion exp Tests" begin
+    s = siteinds("Fermion", 2; conserve_qns=true)
+
+    # Matrix test
+    id_tensor = op("I", s[1])
+    @test id_tensor ≈ exp(0.0 * id_tensor)
+
+    # Tensor test
+    id_tensor = op("I", s[1]) * op("I", s[2])
+    @test id_tensor ≈ exp(0.0 * id_tensor)
+
+    # Permute and test again
+    id_tensor = permute(id_tensor, s[2], s[1], s[2]', s[1]')
+    @test id_tensor ≈ exp(0.0 * id_tensor)
+
+    # Explicitly passing indices in different, valid orders
+    @test id_tensor ≈ exp(0.0 * id_tensor, (s[1]', s[2]'), (dag(s[1]), dag(s[2])))
+    @test id_tensor ≈ exp(0.0 * id_tensor, (s[2]', s[1]'), (dag(s[2]), dag(s[1])))
+    @test id_tensor ≈ exp(0.0 * id_tensor, (dag(s[1]), dag(s[2])), (s[1]', s[2]'))
+    @test id_tensor ≈ exp(0.0 * id_tensor, (dag(s[2]), dag(s[1])), (s[2]', s[1]'))
+
+    # Check wrong index ordering fails (i.e. we are actually paying attention to it)
+    @test norm(id_tensor - exp(0.0 * id_tensor, (dag(s[1]), dag(s[2])), (s[2]', s[1]'))) > 1
+  end
+
   ITensors.disable_auto_fermion()
 end