diff --git a/src/plonk.rs b/src/plonk.rs
index ac5ac13..6def55e 100644
--- a/src/plonk.rs
+++ b/src/plonk.rs
@@ -50,6 +50,7 @@ pub struct Proof<C: CurveAffine> {
     permutation_product_inv_evals: Vec<C::Scalar>,
     permutation_evals: Vec<Vec<C::Scalar>>,
     advice_evals: Vec<C::Scalar>,
+    aux_evals: Vec<C::Scalar>,
     fixed_evals: Vec<C::Scalar>,
     h_evals: Vec<C::Scalar>,
     f_commitment: C,
@@ -91,8 +92,8 @@ fn hash_point<C: CurveAffine, H: Hasher<C::Base>>(
 
 #[test]
 fn test_proving() {
-    use crate::arithmetic::{EqAffine, Field, Fp, Fq};
-    use crate::poly::commitment::Params;
+    use crate::arithmetic::{Curve, EqAffine, Field, Fp, Fq};
+    use crate::poly::commitment::{Blind, Params};
     use crate::transcript::DummyHash;
     use std::marker::PhantomData;
     const K: u32 = 5;
@@ -115,6 +116,7 @@ fn test_proving() {
         sb: FixedWire,
         sc: FixedWire,
         sm: FixedWire,
+        sp: FixedWire,
 
         perm: usize,
         perm2: usize,
@@ -128,6 +130,9 @@ fn test_proving() {
         where
             F: FnOnce() -> Result<(FF, FF, FF), Error>;
         fn copy(&mut self, a: Variable, b: Variable) -> Result<(), Error>;
+        fn public_input<F>(&mut self, f: F) -> Result<Variable, Error>
+        where
+            F: FnOnce() -> Result<FF, Error>;
     }
 
     struct MyCircuit<F: Field> {
@@ -248,6 +253,18 @@ fn test_proving() {
             self.cs
                 .copy(self.config.perm2, left_wire, left.1, right_wire, right.1)
         }
+        fn public_input<F>(&mut self, f: F) -> Result<Variable, Error>
+        where
+            F: FnOnce() -> Result<FF, Error>,
+        {
+            let index = self.current_gate;
+            self.current_gate += 1;
+            self.cs.assign_advice(self.config.a, index, || f())?;
+            self.cs
+                .assign_fixed(self.config.sp, index, || Ok(FF::one()))?;
+
+            Ok(Variable(self.config.a, index))
+        }
     }
 
     impl<F: Field> Circuit<F> for MyCircuit<F> {
@@ -260,6 +277,7 @@ fn test_proving() {
             let sf = meta.fixed_wire();
             let c = meta.advice_wire();
             let d = meta.advice_wire();
+            let p = meta.aux_wire();
 
             let perm = meta.permutation(&[a, b, c]);
             let perm2 = meta.permutation(&[a, b, c]);
@@ -268,6 +286,7 @@ fn test_proving() {
             let sa = meta.fixed_wire();
             let sb = meta.fixed_wire();
             let sc = meta.fixed_wire();
+            let sp = meta.fixed_wire();
 
             meta.create_gate(|meta| {
                 let d = meta.query_advice(d, 1);
@@ -285,6 +304,14 @@ fn test_proving() {
                 a.clone() * sa + b.clone() * sb + a * b * sm + (c * sc * (-F::one())) + sf * (d * e)
             });
 
+            meta.create_gate(|meta| {
+                let a = meta.query_advice(a, 0);
+                let p = meta.query_aux(p, 0);
+                let sp = meta.query_fixed(sp, 0);
+
+                sp * (a + p * (-F::one()))
+            });
+
             PLONKConfig {
                 a,
                 b,
@@ -295,6 +322,7 @@ fn test_proving() {
                 sb,
                 sc,
                 sm,
+                sp,
                 perm,
                 perm2,
             }
@@ -307,6 +335,8 @@ fn test_proving() {
         ) -> Result<(), Error> {
             let mut cs = StandardPLONK::new(cs, config);
 
+            let _ = cs.public_input(|| Ok(F::one() + F::one()))?;
+
             for _ in 0..10 {
                 let mut a_squared = None;
                 let (a0, _, c0) = cs.raw_multiply(|| {
@@ -342,14 +372,26 @@ fn test_proving() {
     // Initialize the SRS
     let srs = SRS::generate(&params, &empty_circuit).expect("SRS generation should not fail");
 
+    let mut pubinputs = srs.domain.empty_lagrange();
+    pubinputs[0] = Fp::one();
+    pubinputs[0] += Fp::one();
+    let pubinput = params
+        .commit_lagrange(&pubinputs, Blind(Field::zero()))
+        .to_affine();
+
     for _ in 0..100 {
         // Create a proof
-        let proof = Proof::create::<DummyHash<Fq>, DummyHash<Fp>, _>(&params, &srs, &circuit)
-            .expect("proof generation should not fail");
+        let proof = Proof::create::<DummyHash<Fq>, DummyHash<Fp>, _>(
+            &params,
+            &srs,
+            &circuit,
+            &[pubinputs.clone()],
+        )
+        .expect("proof generation should not fail");
 
         let msm = params.empty_msm();
         let guard = proof
-            .verify::<DummyHash<Fq>, DummyHash<Fp>>(&params, &srs, msm)
+            .verify::<DummyHash<Fq>, DummyHash<Fp>>(&params, &srs, msm, &[pubinput])
             .unwrap();
         {
             let msm = guard.clone().use_challenges();
@@ -363,7 +405,7 @@ fn test_proving() {
         let msm = guard.clone().use_challenges();
         assert!(msm.clone().is_zero());
         let guard = proof
-            .verify::<DummyHash<Fq>, DummyHash<Fp>>(&params, &srs, msm)
+            .verify::<DummyHash<Fq>, DummyHash<Fp>>(&params, &srs, msm, &[pubinput])
             .unwrap();
         {
             let msm = guard.clone().use_challenges();
diff --git a/src/plonk/circuit.rs b/src/plonk/circuit.rs
index 8ec8e3d..5330a1b 100644
--- a/src/plonk/circuit.rs
+++ b/src/plonk/circuit.rs
@@ -14,6 +14,10 @@ pub struct FixedWire(pub usize);
 #[derive(Copy, Clone, Debug, Eq, PartialEq, Hash)]
 pub struct AdviceWire(pub usize);
 
+/// This represents a wire which has an externally assigned value
+#[derive(Copy, Clone, Debug, Eq, PartialEq, Hash)]
+pub struct AuxWire(pub usize);
+
 /// This trait allows a [`Circuit`] to direct some backend to assign a witness
 /// for a constraint system.
 pub trait Assignment<F: Field> {
@@ -68,6 +72,8 @@ pub enum Expression<F> {
     Fixed(usize),
     /// This is an advice (witness) wire queried at a certain relative location
     Advice(usize),
+    /// This is an auxiliary (external) wire queried at a certain relative location
+    Aux(usize),
     /// This is the sum of two polynomials
     Sum(Box<Expression<F>>, Box<Expression<F>>),
     /// This is the product of two polynomials
@@ -83,6 +89,7 @@ impl<F: Field> Expression<F> {
         &self,
         fixed_wire: &impl Fn(usize) -> T,
         advice_wire: &impl Fn(usize) -> T,
+        aux_wire: &impl Fn(usize) -> T,
         sum: &impl Fn(T, T) -> T,
         product: &impl Fn(T, T) -> T,
         scaled: &impl Fn(T, F) -> T,
@@ -90,18 +97,19 @@ impl<F: Field> Expression<F> {
         match self {
             Expression::Fixed(index) => fixed_wire(*index),
             Expression::Advice(index) => advice_wire(*index),
+            Expression::Aux(index) => aux_wire(*index),
             Expression::Sum(a, b) => {
-                let a = a.evaluate(fixed_wire, advice_wire, sum, product, scaled);
-                let b = b.evaluate(fixed_wire, advice_wire, sum, product, scaled);
+                let a = a.evaluate(fixed_wire, advice_wire, aux_wire, sum, product, scaled);
+                let b = b.evaluate(fixed_wire, advice_wire, aux_wire, sum, product, scaled);
                 sum(a, b)
             }
             Expression::Product(a, b) => {
-                let a = a.evaluate(fixed_wire, advice_wire, sum, product, scaled);
-                let b = b.evaluate(fixed_wire, advice_wire, sum, product, scaled);
+                let a = a.evaluate(fixed_wire, advice_wire, aux_wire, sum, product, scaled);
+                let b = b.evaluate(fixed_wire, advice_wire, aux_wire, sum, product, scaled);
                 product(a, b)
             }
             Expression::Scaled(a, f) => {
-                let a = a.evaluate(fixed_wire, advice_wire, sum, product, scaled);
+                let a = a.evaluate(fixed_wire, advice_wire, aux_wire, sum, product, scaled);
                 scaled(a, *f)
             }
         }
@@ -112,6 +120,7 @@ impl<F: Field> Expression<F> {
         match self {
             Expression::Fixed(_) => 1,
             Expression::Advice(_) => 1,
+            Expression::Aux(_) => 1,
             Expression::Sum(a, b) => max(a.degree(), b.degree()),
             Expression::Product(a, b) => a.degree() + b.degree(),
             Expression::Scaled(poly, _) => poly.degree(),
@@ -151,8 +160,10 @@ pub(crate) struct PointIndex(pub usize);
 pub struct ConstraintSystem<F> {
     pub(crate) num_fixed_wires: usize,
     pub(crate) num_advice_wires: usize,
+    pub(crate) num_aux_wires: usize,
     pub(crate) gates: Vec<Expression<F>>,
     pub(crate) advice_queries: Vec<(AdviceWire, Rotation)>,
+    pub(crate) aux_queries: Vec<(AuxWire, Rotation)>,
     pub(crate) fixed_queries: Vec<(FixedWire, Rotation)>,
 
     // Mapping from a witness vector rotation to the index in the point vector.
@@ -176,9 +187,11 @@ impl<F: Field> Default for ConstraintSystem<F> {
         ConstraintSystem {
             num_fixed_wires: 0,
             num_advice_wires: 0,
+            num_aux_wires: 0,
             gates: vec![],
             fixed_queries: Vec::new(),
             advice_queries: Vec::new(),
+            aux_queries: Vec::new(),
             rotations,
             permutations: Vec::new(),
         }
@@ -255,6 +268,32 @@ impl<F: Field> ConstraintSystem<F> {
         Expression::Advice(self.query_advice_index(wire, at))
     }
 
+    fn query_aux_index(&mut self, wire: AuxWire, at: i32) -> usize {
+        let at = Rotation(at);
+        {
+            let len = self.rotations.len();
+            self.rotations.entry(at).or_insert(PointIndex(len));
+        }
+
+        // Return existing query, if it exists
+        for (index, aux_query) in self.aux_queries.iter().enumerate() {
+            if aux_query == &(wire, at) {
+                return index;
+            }
+        }
+
+        // Make a new query
+        let index = self.aux_queries.len();
+        self.aux_queries.push((wire, at));
+
+        index
+    }
+
+    /// Query an auxiliary wire at a relative position
+    pub fn query_aux(&mut self, wire: AuxWire, at: i32) -> Expression<F> {
+        Expression::Aux(self.query_aux_index(wire, at))
+    }
+
     /// Create a new gate
     pub fn create_gate(&mut self, f: impl FnOnce(&mut Self) -> Expression<F>) {
         let poly = f(self);
@@ -274,4 +313,11 @@ impl<F: Field> ConstraintSystem<F> {
         self.num_advice_wires += 1;
         tmp
     }
+
+    /// Allocate a new auxiliary wire
+    pub fn aux_wire(&mut self) -> AuxWire {
+        let tmp = AuxWire(self.num_aux_wires);
+        self.num_aux_wires += 1;
+        tmp
+    }
 }
diff --git a/src/plonk/prover.rs b/src/plonk/prover.rs
index 93c6942..7c7d514 100644
--- a/src/plonk/prover.rs
+++ b/src/plonk/prover.rs
@@ -24,7 +24,12 @@ impl<C: CurveAffine> Proof<C> {
         params: &Params<C>,
         srs: &SRS<C>,
         circuit: &ConcreteCircuit,
+        aux: &[Polynomial<C::Scalar, LagrangeCoeff>],
     ) -> Result<Self, Error> {
+        if aux.len() != srs.cs.num_aux_wires {
+            return Err(Error::IncompatibleParams);
+        }
+
         struct WitnessCollection<F: Field> {
             advice: Vec<Polynomial<F, LagrangeCoeff>>,
             _marker: std::marker::PhantomData<F>,
@@ -88,6 +93,38 @@ impl<C: CurveAffine> Proof<C> {
         // Create a transcript for obtaining Fiat-Shamir challenges.
         let mut transcript = HBase::init(C::Base::one());
 
+        // Compute commitments to aux wire polynomials
+        let aux_commitments_projective: Vec<_> = aux
+            .iter()
+            .map(|poly| params.commit_lagrange(poly, Blind(C::Scalar::zero()))) // TODO: bad blind?
+            .collect();
+        let mut aux_commitments = vec![C::zero(); aux_commitments_projective.len()];
+        C::Projective::batch_to_affine(&aux_commitments_projective, &mut aux_commitments);
+        let aux_commitments = aux_commitments;
+        drop(aux_commitments_projective);
+
+        for commitment in &aux_commitments {
+            hash_point(&mut transcript, commitment)?;
+        }
+
+        let aux_polys: Vec<_> = aux
+            .clone()
+            .into_iter()
+            .map(|poly| {
+                let lagrange_vec = domain.lagrange_from_vec(poly.to_vec());
+                domain.lagrange_to_coeff(lagrange_vec)
+            })
+            .collect();
+
+        let aux_cosets: Vec<_> = meta
+            .aux_queries
+            .iter()
+            .map(|&(wire, at)| {
+                let poly = aux_polys[wire.0].clone();
+                domain.coeff_to_extended(poly, at)
+            })
+            .collect();
+
         // Compute commitments to advice wire polynomials
         let advice_blinds: Vec<_> = witness
             .advice
@@ -253,6 +290,7 @@ impl<C: CurveAffine> Proof<C> {
             let evaluation = poly.evaluate(
                 &|index| srs.fixed_cosets[index].clone(),
                 &|index| advice_cosets[index].clone(),
+                &|index| aux_cosets[index].clone(),
                 &|a, b| a + &b,
                 &|a, b| a * &b,
                 &|a, scalar| a * scalar,
@@ -355,6 +393,12 @@ impl<C: CurveAffine> Proof<C> {
             .map(|&(wire, at)| eval_polynomial(&advice_polys[wire.0], domain.rotate_omega(x_3, at)))
             .collect();
 
+        let aux_evals: Vec<_> = meta
+            .aux_queries
+            .iter()
+            .map(|&(wire, at)| eval_polynomial(&aux_polys[wire.0], domain.rotate_omega(x_3, at)))
+            .collect();
+
         let fixed_evals: Vec<_> = meta
             .fixed_queries
             .iter()
@@ -396,6 +440,7 @@ impl<C: CurveAffine> Proof<C> {
         // Hash each advice evaluation
         for eval in advice_evals
             .iter()
+            .chain(aux_evals.iter())
             .chain(fixed_evals.iter())
             .chain(h_evals.iter())
             .chain(permutation_product_evals.iter())
@@ -451,6 +496,17 @@ impl<C: CurveAffine> Proof<C> {
                 );
             }
 
+            for (query_index, &(wire, ref at)) in meta.aux_queries.iter().enumerate() {
+                let point_index = (*meta.rotations.get(at).unwrap()).0;
+
+                accumulate(
+                    point_index,
+                    &aux_polys[wire.0],
+                    Blind(C::Scalar::zero()),
+                    aux_evals[query_index],
+                );
+            }
+
             for (query_index, &(wire, ref at)) in meta.fixed_queries.iter().enumerate() {
                 let point_index = (*meta.rotations.get(at).unwrap()).0;
 
@@ -595,6 +651,7 @@ impl<C: CurveAffine> Proof<C> {
             permutation_evals,
             advice_evals,
             fixed_evals,
+            aux_evals,
             h_evals,
             f_commitment,
             q_evals,
diff --git a/src/plonk/verifier.rs b/src/plonk/verifier.rs
index a9ae739..0200e93 100644
--- a/src/plonk/verifier.rs
+++ b/src/plonk/verifier.rs
@@ -13,7 +13,16 @@ impl<'a, C: CurveAffine> Proof<C> {
         params: &'a Params<C>,
         srs: &SRS<C>,
         mut msm: MSM<'a, C>,
+        aux_commitments: &[C],
     ) -> Result<Guard<'a, C>, Error> {
+        // Check that aux_commitments matches the expected number of aux_wires
+        // and self.aux_evals
+        if aux_commitments.len() != srs.cs.num_aux_wires
+            || self.aux_evals.len() != srs.cs.num_aux_wires
+        {
+            return Err(Error::IncompatibleParams);
+        }
+
         // Scale the MSM by a random factor to ensure that if the existing MSM
         // has is_zero() == false then this argument won't be able to interfere
         // with it to make it true, with high probability.
@@ -22,10 +31,14 @@ impl<'a, C: CurveAffine> Proof<C> {
         // Create a transcript for obtaining Fiat-Shamir challenges.
         let mut transcript = HBase::init(C::Base::one());
 
+        // Hash the aux (external) commitments into the transcript
+        for commitment in aux_commitments {
+            hash_point(&mut transcript, commitment)?;
+        }
+
         // Hash the prover's advice commitments into the transcript
         for commitment in &self.advice_commitments {
-            hash_point(&mut transcript, commitment)
-                .expect("proof cannot contain points at infinity");
+            hash_point(&mut transcript, commitment)?;
         }
 
         // Sample x_0 challenge
@@ -36,7 +49,7 @@ impl<'a, C: CurveAffine> Proof<C> {
 
         // Hash each permutation product commitment
         for c in &self.permutation_product_commitments {
-            hash_point(&mut transcript, c).expect("proof cannot contain points at infinity");
+            hash_point(&mut transcript, c)?;
         }
 
         // Sample x_2 challenge, which keeps the gates linearly independent.
@@ -44,7 +57,7 @@ impl<'a, C: CurveAffine> Proof<C> {
 
         // Obtain a commitment to h(X) in the form of multiple pieces of degree n - 1
         for c in &self.h_commitments {
-            hash_point(&mut transcript, c).expect("proof cannot contain points at infinity");
+            hash_point(&mut transcript, c)?;
         }
 
         // Sample x_3 challenge, which is used to ensure the circuit is
@@ -59,6 +72,7 @@ impl<'a, C: CurveAffine> Proof<C> {
         for eval in self
             .advice_evals
             .iter()
+            .chain(self.aux_evals.iter())
             .chain(self.fixed_evals.iter())
             .chain(self.h_evals.iter())
             .chain(self.permutation_product_evals.iter())
@@ -80,6 +94,7 @@ impl<'a, C: CurveAffine> Proof<C> {
             let evaluation: C::Scalar = poly.evaluate(
                 &|index| self.fixed_evals[index],
                 &|index| self.advice_evals[index],
+                &|index| self.aux_evals[index],
                 &|a, b| a + &b,
                 &|a, b| a * &b,
                 &|a, scalar| a * &scalar,
@@ -172,6 +187,15 @@ impl<'a, C: CurveAffine> Proof<C> {
                 );
             }
 
+            for (query_index, &(wire, ref at)) in srs.cs.aux_queries.iter().enumerate() {
+                let point_index = (*srs.cs.rotations.get(at).unwrap()).0;
+                accumulate(
+                    point_index,
+                    aux_commitments[wire.0],
+                    self.aux_evals[query_index],
+                );
+            }
+
             for (query_index, &(wire, ref at)) in srs.cs.fixed_queries.iter().enumerate() {
                 let point_index = (*srs.cs.rotations.get(at).unwrap()).0;
                 accumulate(
@@ -222,8 +246,7 @@ impl<'a, C: CurveAffine> Proof<C> {
         let x_5: C::Scalar = get_challenge_scalar(Challenge(transcript.squeeze().get_lower_128()));
 
         // Obtain the commitment to the multi-point quotient polynomial f(X).
-        hash_point(&mut transcript, &self.f_commitment)
-            .expect("proof cannot contain points at infinity");
+        hash_point(&mut transcript, &self.f_commitment)?;
 
         // Sample a challenge x_6 for checking that f(X) was committed to
         // correctly.