diff --git a/Project.toml b/Project.toml
index cee56bc..dcb64cd 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "HarmonicOrthogonalPolynomials"
 uuid = "e416a80e-9640-42f3-8df8-80a93ca01ea5"
 authors = ["Sheehan Olver <solver@mac.com>"]
-version = "0.6.3"
+version = "0.7"
 
 [deps]
 BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
@@ -20,13 +20,13 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 [compat]
 BlockArrays = "1.0"
 BlockBandedMatrices = "0.13"
-ClassicalOrthogonalPolynomials = "0.13, 0.14"
-ContinuumArrays = "0.18"
+ClassicalOrthogonalPolynomials = "0.15"
+ContinuumArrays = "0.19"
 DomainSets = "0.7"
 FastTransforms = "0.15, 0.16, 0.17"
-InfiniteArrays = "0.14, 0.15"
+InfiniteArrays = "0.15"
 IntervalSets = "0.7"
-QuasiArrays = "0.11"
+QuasiArrays = "0.12"
 SpecialFunctions = "1, 2"
 StaticArrays = "1"
 julia = "1.10"
diff --git a/src/HarmonicOrthogonalPolynomials.jl b/src/HarmonicOrthogonalPolynomials.jl
index 68a79e4..7af6175 100644
--- a/src/HarmonicOrthogonalPolynomials.jl
+++ b/src/HarmonicOrthogonalPolynomials.jl
@@ -5,15 +5,17 @@ import Base: OneTo, oneto, axes, getindex, convert, to_indices, tail, eltype, *,
 import BlockArrays: block, blockindex, unblock, BlockSlice
 import DomainSets: indomain, Sphere
 import LinearAlgebra: norm, factorize
-import QuasiArrays: to_quasi_index, SubQuasiArray, *
-import ContinuumArrays: TransformFactorization, @simplify, ProjectionFactorization, plan_grid_transform, plan_transform, grid, grid_layout, plotgrid_layout, AbstractBasisLayout, MemoryLayout
+import QuasiArrays: to_quasi_index, SubQuasiArray, *, AbstractQuasiVecOrMat
+import ContinuumArrays: TransformFactorization, @simplify, ProjectionFactorization, plan_grid_transform, plan_transform, grid, grid_layout, plotgrid_layout,
+                        AbstractBasisLayout, MemoryLayout, abslaplacian, laplacian, AbstractDifferentialQuasiMatrix, operatorcall, similaroperator, SubBasisLayout,
+                        ApplyLayout, arguments, ExpansionLayout
 import ClassicalOrthogonalPolynomials: checkpoints, _sum, cardinality, increasingtruncations
 import BlockBandedMatrices: BlockRange1, _BandedBlockBandedMatrix
 import FastTransforms: Plan, interlace
 import QuasiArrays: LazyQuasiMatrix, LazyQuasiArrayStyle
 import InfiniteArrays: InfStepRange, RangeCumsum
 
-export SphericalHarmonic, UnitSphere, SphericalCoordinate, RadialCoordinate, Block, associatedlegendre, RealSphericalHarmonic, sphericalharmonicy, Laplacian, AbsLaplacianPower, abs, -, ^, AngularMomentum
+export SphericalHarmonic, UnitSphere, SphericalCoordinate, RadialCoordinate, Block, associatedlegendre, RealSphericalHarmonic, sphericalharmonicy, abs, -, ^, AngularMomentum, Laplacian, AbsLaplacian
 cardinality(::Sphere) = ℵ₁
 
 include("multivariateops.jl")
diff --git a/src/angularmomentum.jl b/src/angularmomentum.jl
index b2faa59..6a67897 100644
--- a/src/angularmomentum.jl
+++ b/src/angularmomentum.jl
@@ -1,25 +1,54 @@
-#########
-# AngularMomentum
-# Applies the partial derivative with respect to the last angular variable in the coordinate system.
-# For example, in polar coordinates (r, θ) in ℝ² or cylindrical coordinates (r, θ, z) in ℝ³, we apply ∂ / ∂θ = (x ∂ / ∂y - y ∂ / ∂x).
-# In spherical coordinates (ρ, θ, φ) in ℝ³, we apply ∂ / ∂φ = (x ∂ / ∂y - y ∂ / ∂x).
-#########
-
-struct AngularMomentum{T,Ax<:Inclusion} <: LazyQuasiMatrix{T}
+"""
+    AngularMomentum
+    
+Applies the partial derivative with respect to the last angular variable in the coordinate system.
+For example, in polar coordinates (r, θ) in ℝ² or cylindrical coordinates (r, θ, z) in ℝ³, we apply ∂ / ∂θ = (x ∂ / ∂y - y ∂ / ∂x).
+In spherical coordinates (ρ, θ, φ) in ℝ³, we apply ∂ / ∂φ = (x ∂ / ∂y - y ∂ / ∂x).
+"""
+struct AngularMomentum{T,Ax<:Inclusion,Order} <: AbstractDifferentialQuasiMatrix{T}
     axis::Ax
+    order::Order
 end
 
-AngularMomentum{T}(axis::Inclusion) where T = AngularMomentum{T,typeof(axis)}(axis)
-AngularMomentum{T}(domain) where T = AngularMomentum{T}(Inclusion(domain))
-AngularMomentum(axis) = AngularMomentum{eltype(eltype(axis))}(axis)
+AngularMomentum{T, D}(axis::D, order) where {T,D<:Inclusion} = AngularMomentum{T,D,typeof(order)}(axis, order)
+AngularMomentum{T, D}(axis::D) where {T,D<:Inclusion} = AngularMomentum{T,D,Nothing}(axis, nothing)
+AngularMomentum{T}(axis::Inclusion, order...) where T = AngularMomentum{T,typeof(axis)}(axis, order...)
+AngularMomentum{T}(domain, order...) where T = AngularMomentum{T}(Inclusion(domain), order...)
+AngularMomentum(domain, order...) = AngularMomentum(Inclusion(domain), order...)
+AngularMomentum(ax::AbstractQuasiVector{T}, order...) where T = AngularMomentum{eltype(eltype(ax))}(ax, order...)
+AngularMomentum(L::AbstractQuasiMatrix, order...) = AngularMomentum(axes(L,1), order...)
 
-axes(A::AngularMomentum) = (A.axis, A.axis)
-==(a::AngularMomentum, b::AngularMomentum) = a.axis == b.axis
-copy(A::AngularMomentum) = AngularMomentum(copy(A.axis))
 
-^(A::AngularMomentum, k::Integer) = ApplyQuasiArray(^, A, k)
+operatorcall(::AngularMomentum) = angularmomentum
+similaroperator(D::AngularMomentum, ax, ord) = AngularMomentum{eltype(D)}(ax, ord)
 
-@simplify function *(A::AngularMomentum, P::SphericalHarmonic)
+simplifiable(::typeof(*), A::AngularMomentum, B::AngularMomentum) = Val(true)
+*(D1::AngularMomentum, D2::AngularMomentum) = similaroperator(convert(AbstractQuasiMatrix{promote_type(eltype(D1),eltype(D2))}, D1), D1.axis, operatororder(D1)+operatororder(D2))
+
+angularmomentum(A, order...; dims...) = angularmomentum_layout(MemoryLayout(A), A, order...; dims...)
+
+angularmomentum_layout(::AbstractBasisLayout, Vm, order...; dims...) = error("Overload angularmomentum(::$(typeof(Vm)))")
+function angularmomentum_layout(::AbstractBasisLayout, a, order::Int; dims...)
+    order < 0 && throw(ArgumentError("order must be non-negative"))
+    order == 0 && return a
+    isone(order) ? angularmomentum(a) : angularmomentum(angularmomentum(a), order-1)
+end
+
+angularmomentum_layout(::ExpansionLayout, A, order...; dims...) = angularmomentum_layout(ApplyLayout{typeof(*)}(), A, order...; dims...)
+
+
+function angularmomentum_layout(::SubBasisLayout, Vm, order...; dims::Integer=1)
+    dims == 1 || error("not implemented")
+    angularmomentum(parent(Vm), order...)[:,parentindices(Vm)[2]]
+end
+
+function angularmomentum_layout(LAY::ApplyLayout{typeof(*)}, V::AbstractQuasiVecOrMat, order...; dims=1)
+    a = arguments(LAY, V)
+    dims == 1 || throw(ArgumentError("cannot take angularmomentum a vector along dimension $dims"))
+    *(angularmomentum(a[1], order...), tail(a)...)
+end
+
+function angularmomentum(P::SphericalHarmonic)
     # Spherical harmonics are the eigenfunctions of the angular momentum operator on the unit sphere
     T = real(eltype(P))
     k = mortar(Base.OneTo.(1:2:∞))
diff --git a/src/laplace.jl b/src/laplace.jl
index 2a6ecb3..345c439 100644
--- a/src/laplace.jl
+++ b/src/laplace.jl
@@ -1,36 +1,20 @@
-# Laplacian
-
-struct Laplacian{T,D} <: LazyQuasiMatrix{T}
-    axis::Inclusion{T,D}
-end
-
-Laplacian{T}(axis::Inclusion{<:Any,D}) where {T,D} = Laplacian{T,D}(axis)
-# Laplacian{T}(domain) where T = Laplacian{T}(Inclusion(domain))
-# Laplacian(axis) = Laplacian{eltype(axis)}(axis)
-
-axes(D::Laplacian) = (D.axis, D.axis)
-==(a::Laplacian, b::Laplacian) = a.axis == b.axis
-copy(D::Laplacian) = Laplacian(copy(D.axis))
-
-@simplify function *(Δ::Laplacian, P::AbstractSphericalHarmonic)
+function abslaplacian(P::AbstractSphericalHarmonic; dims...)
      # Spherical harmonics are the eigenfunctions of the Laplace operator on the unit sphere
-     P * Diagonal(mortar(Fill.((-(0:∞)-(0:∞).^2), 1:2:∞)))
+     P * Diagonal(mortar(Fill.((0:∞)+(0:∞).^2, 1:2:∞)))
 end
 
-# Negative fractional Laplacian (-Δ)^α or equiv. abs(Δ)^α
-
-struct AbsLaplacianPower{T,D,A} <: LazyQuasiMatrix{T}
-    axis::Inclusion{T,D}
-    α::A
+function abslaplacian(P::AbstractSphericalHarmonic, order; dims...)
+    # Spherical harmonics are the eigenfunctions of the Laplace operator on the unit sphere
+    P * Diagonal(mortar(Fill.(((0:∞)+(0:∞).^2) .^ order, 1:2:∞)))
 end
 
-AbsLaplacianPower{T}(axis::Inclusion{<:Any,D},α) where {T,D} = AbsLaplacianPower{T,D,typeof(α)}(axis,α)
-
-axes(D:: AbsLaplacianPower) = (D.axis, D.axis)
-==(a:: AbsLaplacianPower, b:: AbsLaplacianPower) = a.axis == b.axis && a.α == b.α
-copy(D:: AbsLaplacianPower) = AbsLaplacianPower(copy(D.axis), D.α)
 
-abs(Δ::Laplacian) = AbsLaplacianPower(axes(Δ,1),1)
--(Δ::Laplacian) = abs(Δ)
+function laplacian(P::AbstractSphericalHarmonic; dims...)
+    # Spherical harmonics are the eigenfunctions of the Laplace operator on the unit sphere
+    P * Diagonal(mortar(Fill.(-(0:∞)-(0:∞).^2, 1:2:∞)))
+end
 
-^(D::AbsLaplacianPower, k) = AbsLaplacianPower(D.axis, D.α*k)
\ No newline at end of file
+function laplacian(P::AbstractSphericalHarmonic, order::Int; dims...)
+    # Spherical harmonics are the eigenfunctions of the Laplace operator on the unit sphere
+    P * Diagonal(mortar(Fill.((-(0:∞)-(0:∞).^2) .^ order, 1:2:∞)))
+end
diff --git a/src/multivariateops.jl b/src/multivariateops.jl
index ae8aa3e..3e69b4d 100644
--- a/src/multivariateops.jl
+++ b/src/multivariateops.jl
@@ -1,24 +1,3 @@
-
-#########
-# PartialDerivative{k}
-# takes a partial derivative in the k-th index.
-#########
-
-
-struct PartialDerivative{k,T,Ax<:Inclusion} <: LazyQuasiMatrix{T}
-    axis::Ax
-end
-
-PartialDerivative{k,T}(axis::Inclusion) where {k,T} = PartialDerivative{k,T,typeof(axis)}(axis)
-PartialDerivative{k,T}(domain) where {k,T} = PartialDerivative{k,T}(Inclusion(domain))
-PartialDerivative{k}(axis) where k = PartialDerivative{k,eltype(eltype(axis))}(axis)
-
-axes(D::PartialDerivative) = (D.axis, D.axis)
-==(a::PartialDerivative{k}, b::PartialDerivative{k}) where k = a.axis == b.axis
-copy(D::PartialDerivative{k}) where k = PartialDerivative{k}(copy(D.axis))
-
-^(D::PartialDerivative, k::Integer) = ApplyQuasiArray(^, D, k)
-
 abstract type MultivariateOrthogonalPolynomial{d,T} <: Basis{T} end
 const BivariateOrthogonalPolynomial{T} = MultivariateOrthogonalPolynomial{2,T}
 
diff --git a/test/runtests.jl b/test/runtests.jl
index 848feb2..0cd999a 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -229,7 +229,7 @@ end
     RΔ = Laplacian(Rxyz)
     @test SΔ isa Laplacian
     @test RΔ isa Laplacian
-    @test *(SΔ,S) isa ApplyQuasiArray
+    @test SΔ*S isa ApplyQuasiArray
     @test *(RΔ,R) isa ApplyQuasiArray
     @test copy(SΔ) == SΔ == RΔ == copy(RΔ)
     @test axes(SΔ) == axes(RΔ) == (axes(S,1),axes(S,1)) == (axes(R,1),axes(R,1))
@@ -323,8 +323,8 @@ end
     S = SphericalHarmonic()
     xyz = axes(S,1)
     @test Laplacian(xyz) isa Laplacian
-    @test Laplacian(xyz)^2 isa QuasiArrays.ApplyQuasiArray
-    @test Laplacian(xyz)^3 isa QuasiArrays.ApplyQuasiArray
+    @test Laplacian(xyz)^2 isa Laplacian
+    @test Laplacian(xyz)^3 isa Laplacian
     f1  = c -> cos(c.θ)^2
     Δ_f1 = c -> -1-3*cos(2*c.θ)
     Δ2_f1 = c -> 6+18*cos(2*c.θ)
@@ -342,8 +342,8 @@ end
     S = SphericalHarmonic()[:,Block.(Base.OneTo(10))]
     xyz = axes(S,1)
     @test Laplacian(xyz) isa Laplacian
-    @test Laplacian(xyz)^2 isa QuasiArrays.ApplyQuasiArray
-    @test Laplacian(xyz)^3 isa QuasiArrays.ApplyQuasiArray
+    @test Laplacian(xyz)^2 isa Laplacian
+    @test Laplacian(xyz)^3 isa Laplacian
     f1  = c -> cos(c.θ)^2
     Δ_f1 = c -> -1-3*cos(2*c.θ)
     Δ2_f1 = c -> 6+18*cos(2*c.θ)
@@ -361,8 +361,8 @@ end
     S = RealSphericalHarmonic()[:,Block.(Base.OneTo(10))]
     xyz = axes(S,1)
     @test Laplacian(xyz) isa Laplacian
-    @test Laplacian(xyz)^2 isa QuasiArrays.ApplyQuasiArray
-    @test Laplacian(xyz)^3 isa QuasiArrays.ApplyQuasiArray
+    @test Laplacian(xyz)^2 isa Laplacian
+    @test Laplacian(xyz)^3 isa Laplacian
     f1  = c -> cos(c.θ)^2
     Δ_f1 = c -> -1-3*cos(2*c.θ)
     Δ2_f1 = c -> 6+18*cos(2*c.θ)
@@ -381,25 +381,25 @@ end
     α = 1/3
     S = SphericalHarmonic()
     Sxyz = axes(S,1)
-    SΔα = AbsLaplacianPower(Sxyz,α)
+    SΔα = AbsLaplacian(Sxyz,α)
     Δ = Laplacian(Sxyz)
     @test copy(SΔα) == SΔα
-    @test SΔα isa AbsLaplacianPower
+    @test SΔα isa AbsLaplacian
     @test SΔα isa QuasiArrays.LazyQuasiMatrix
     @test axes(SΔα) == (axes(S,1),axes(S,1))
-    @test abs(Δ) == -Δ == AbsLaplacianPower(axes(Δ,1),1)
+    @test abs(Δ) == -Δ == AbsLaplacian(axes(Δ,1),1)
     @test abs(Δ)^α == SΔα
     # Set 2
     α = 7/13
     S = SphericalHarmonic()
     Sxyz = axes(S,1)
-    SΔα = AbsLaplacianPower(Sxyz,α)
+    SΔα = AbsLaplacian(Sxyz,α)
     Δ = Laplacian(Sxyz)
     @test copy(SΔα) == SΔα
-    @test SΔα isa AbsLaplacianPower
+    @test SΔα isa AbsLaplacian
     @test SΔα isa QuasiArrays.LazyQuasiMatrix
     @test axes(SΔα) == (axes(S,1),axes(S,1))
-    @test abs(Δ) == -Δ == AbsLaplacianPower(axes(Δ,1),1)
+    @test abs(Δ) == -Δ == AbsLaplacian(axes(Δ,1),1)
     @test abs(Δ)^α == SΔα
 end
 
@@ -417,9 +417,9 @@ end
 @testset "Angular momentum" begin
     S = SphericalHarmonic()
     R = RealSphericalHarmonic()
-    ∂θ = AngularMomentum(axes(S, 1))
+    ∂θ = AngularMomentum(S)
     @test axes(∂θ) == (axes(S, 1), axes(S, 1))
-    @test ∂θ == AngularMomentum(axes(R, 1)) == AngularMomentum(axes(S, 1).domain)
+    @test ∂θ == AngularMomentum(R) == AngularMomentum(axes(S, 1).domain)
     @test copy(∂θ) ≡ ∂θ
     A = S \ (∂θ * S)
     A2 = S \ (∂θ^2 * S)