diff --git a/docs/src/man/models.md b/docs/src/man/models.md
index 8693d4f..7bdba69 100644
--- a/docs/src/man/models.md
+++ b/docs/src/man/models.md
@@ -15,6 +15,7 @@ heisenberg_XXZ
 heisenberg_XYZ
 bilinear_biquadratic_model
 tj_model
+ashkin_teller
 hubbard_model
 bose_hubbard_model
 quantum_chemistry_hamiltonian
diff --git a/src/MPSKitModels.jl b/src/MPSKitModels.jl
index 64e8064..c6aa42f 100644
--- a/src/MPSKitModels.jl
+++ b/src/MPSKitModels.jl
@@ -45,6 +45,7 @@ export heisenberg_XXX, heisenberg_XXZ, heisenberg_XYZ
 export bilinear_biquadratic_model
 export hubbard_model, bose_hubbard_model
 export tj_model
+export ashkin_teller
 export quantum_chemistry_hamiltonian
 
 export classical_ising
diff --git a/src/models/hamiltonians.jl b/src/models/hamiltonians.jl
index 7d36c28..6d448c9 100644
--- a/src/models/hamiltonians.jl
+++ b/src/models/hamiltonians.jl
@@ -461,4 +461,75 @@ function tj_model(
     end
 end
 
+#===========================================================================================
+    Quantum Ashkin-Teller model
+===========================================================================================#
+"""
+    ashkin_teller([T::Type{<:Number} = ComplexF64], [S = Trivial],
+                  [lattice::AbstractLattice = InfiniteChain(1)];
+                  h = 1.0, J = 1.0, λ = 1.0)
+
+MPO for the hamiltonian of the quantum Ashkin-Teller model. The model is
+defined on a chain of two qubits per site. Writing Pauli operators on each of these qubits
+as ``\\sigma `` and ``\\tau ``, the Hamiltonian reads:
+```math
+H = -h \\sum_i\\bigg(\\sigma_i^x + \\tau_i^x + \\lambda \\sigma_i^x \\tau_i^x \\bigg)
+-J \\sum_{\\langle i,j \\rangle} \\bigg( \\sigma_i^z \\sigma_j^z + \\tau_i^z \\tau_j^z
++\\lambda \\sigma_i^z \\sigma_j^z \\tau_i^z \\tau_j^z \\bigg).
+
+```
+"""
+function ashkin_teller end
+function ashkin_teller(lattice::AbstractLattice; kwargs...)
+    return ashkin_teller(ComplexF64, Trivial, lattice; kwargs...)
+end
+function ashkin_teller(symmetry::Type{<:Sector}; kwargs...)
+    return ashkin_teller(ComplexF64, symmetry; kwargs...)
+end
+function ashkin_teller(T::Type{<:Number}, lattice::AbstractLattice; kwargs...)
+    return ashkin_teller(T, Trivial, lattice; kwargs...)
+end
+function ashkin_teller(
+        T::Type{<:Number} = ComplexF64, S::Type{<:Sector} = Trivial,
+        lattice::AbstractLattice = InfiniteChain(1);
+        h = 1.0, J = 1.0, λ = 1.0
+    )
+    S₁, S₂ = _ashin_teller_decompose_symmetry(S)
+
+    # component tensors
+    X₁ = σˣ(T, S₁)
+    Z₁Z₁ = σᶻᶻ(T, S₁)
+    I₁ = id(T, domain(X₁))
+    I₁I₁ = id(T, domain(Z₁Z₁))
+
+    Z₂Z₂ = σᶻᶻ(T, S₂)
+    X₂ = σˣ(T, S₂)
+    I₂ = id(T, domain(X₂))
+    I₂I₂ = id(T, domain(Z₂Z₂))
+
+    # Single site operators
+    XI = X₁ ⊠ I₂
+    IX = I₁ ⊠ X₂
+    XX = X₁ ⊠ X₂
+    F = isometry(Int, fuse(domain(XX)) ← domain(XX))
+    onsite = F * (-h * (XI + IX + λ * XX)) * F'
+
+    # Nearest-neighbour terms
+    @tensor FF[-1 -2; -3 -5 -4 -6] := F[-1; -3 -4] * F[-2; -5 -6] # note permutation!
+    ZIZI = FF * (Z₁Z₁ ⊠ I₂I₂) * FF'
+    IZIZ = FF * (I₁I₁ ⊠ Z₂Z₂) * FF'
+    ZZZZ = FF * (Z₁Z₁ ⊠ Z₂Z₂) * FF'
+    twosite = -J * (ZIZI + IZIZ + λ * ZZZZ)
+
+    return @mpoham sum(vertices(lattice)) do i
+        return onsite{i}
+    end + sum(nearest_neighbours(lattice)) do (i, j)
+        return twosite{i, j}
+    end
+end
+
+_ashin_teller_decompose_symmetry(::Type{Trivial}) = (Trivial, Trivial)
+_ashin_teller_decompose_symmetry(::Type{ProductSector{Tuple{A, B}}}) where {A <: Union{Trivial, Z2Irrep}, B <: Union{Trivial, Z2Irrep}} = (A, B)
+_ashin_teller_decompose_symmetry(T) = error("Ashkin-Teller model not implemented for symmetry $T")
+
 # TODO: add (hardcore) bosonic t-J model (https://arxiv.org/abs/2409.15424)
diff --git a/src/utility.jl b/src/utility.jl
index a80fdfa..ad893df 100644
--- a/src/utility.jl
+++ b/src/utility.jl
@@ -57,3 +57,15 @@ function split_twosite(O::AbstractTensorMap{<:Any, <:Any, 2, 2})
     @plansor R[a p'; p] := sqrtS[a; 1] * V[1; p p']
     return L, R
 end
+
+#===========================================================================================
+    Other utility
+===========================================================================================#
+
+"""
+    flip_charge(charge, flip)
+
+    Take the product of charge with flip.
+    Works only for abelian symmetries.
+"""
+flip_charge(charge, flip) = only(charge ⊗ sectortype(charge)(flip))
diff --git a/test/ashkin_teller.jl b/test/ashkin_teller.jl
new file mode 100644
index 0000000..e0ca98f
--- /dev/null
+++ b/test/ashkin_teller.jl
@@ -0,0 +1,36 @@
+using Test
+using TensorKit
+using MPSKit
+using MPSKitModels
+
+# Setup
+
+gammas = [
+    pi / 2, pi / 3, pi / 4, pi / 6, 0,
+]
+# Exact GS energy density via mapping to spin 1/2 XXZ and Bethe ansatz
+# https://www.sciencedirect.com/science/article/pii/0003491688900152
+E0s = [
+    -8 / pi, -12 / 4, -4 * (sqrt(2) / pi + 1 / (2 * sqrt(2))),
+    -4 * (1 / pi + 11 / (12 * sqrt(3))), 2 - 8 * log(2),
+]
+alg = VUMPS(; maxiter = 100, verbosity = 0)
+
+# Test
+
+@testset "Ashkin-Teller" for (gamma, E0) in zip(gammas, E0s)
+    H = ashkin_teller(h = 1, J = 1, λ = cos(gamma))
+    Ψ = InfiniteMPS(physicalspace(H), [ComplexSpace(24)])
+    @test imag(expectation_value(Ψ, H)) ≈ 0 atol = 1.0e-12
+    Ψ0, _ = find_groundstate(Ψ, H, alg)
+    @test real(expectation_value(Ψ0, H)) ≈ E0 atol = 1.0e-3
+end
+
+@testset "Ashkin-Teller (Z2 ⊠ Z2)" for (gamma, E0) in zip(gammas, E0s)
+    H = ashkin_teller(Z2Irrep ⊠ Z2Irrep; h = 1, J = 1, λ = cos(gamma))
+    V = spacetype(H)(c => 6 for c in values(sectortype(H)))
+    Ψ = InfiniteMPS(physicalspace(H), [V])
+    @test imag(expectation_value(Ψ, H)) ≈ 0 atol = 1.0e-12
+    Ψ0, _ = find_groundstate(Ψ, H, alg)
+    @test real(expectation_value(Ψ0, H)) ≈ E0 atol = 1.0e-3
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 9fc5e38..41a5576 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -45,6 +45,10 @@ end
     include("quantum_potts.jl")
 end
 
+@testset "quantum Ashkin-Teller model" begin
+    include("ashkin_teller.jl")
+end
+
 @testset "classical ising model" begin
     include("classical_ising.jl")
 end