diff --git a/qiskit/aqua/operators/evolutions/evolved_op.py b/qiskit/aqua/operators/evolutions/evolved_op.py
index 2dbb2340fe84..e6fe00991866 100644
--- a/qiskit/aqua/operators/evolutions/evolved_op.py
+++ b/qiskit/aqua/operators/evolutions/evolved_op.py
@@ -86,13 +86,23 @@ def tensor(self, other: OperatorBase) -> OperatorBase:
 
         return TensoredOp([self, other])
 
-    def compose(self, other: OperatorBase) -> OperatorBase:
-        other = self._check_zero_for_composition_and_expand(other)
+    def _expand_dim(self, num_qubits: int) -> OperatorBase:
+        from qiskit.aqua.operators import I
+        return self.tensor(I ^ num_qubits)
 
+    def permute(self, permutation: List[int]) -> OperatorBase:
+        return EvolvedOp(self.primitive.permute(permutation), coeff=self.coeff)  # type: ignore
+
+    def compose(self, other: OperatorBase,
+                permutation: Optional[List[int]] = None, front: bool = False) -> OperatorBase:
+
+        new_self, other = self._expand_shorter_operator_and_permute(other, permutation)
+        if front:
+            return other.compose(new_self)
         if isinstance(other, ComposedOp):
-            return ComposedOp([self] + other.oplist)  # type: ignore
+            return ComposedOp([new_self] + other.oplist)
 
-        return ComposedOp([self, other])
+        return ComposedOp([new_self, other])
 
     def __str__(self) -> str:
         prim_str = str(self.primitive)
diff --git a/qiskit/aqua/operators/list_ops/composed_op.py b/qiskit/aqua/operators/list_ops/composed_op.py
index 4e41fbe51c96..fb6419c01aa3 100644
--- a/qiskit/aqua/operators/list_ops/composed_op.py
+++ b/qiskit/aqua/operators/list_ops/composed_op.py
@@ -12,7 +12,7 @@
 
 """ ComposedOp Class """
 
-from typing import List, Union, cast
+from typing import List, Union, cast, Optional
 from functools import reduce, partial
 import numpy as np
 
@@ -82,24 +82,31 @@ def to_circuit(self) -> QuantumCircuit:
     def adjoint(self) -> OperatorBase:
         return ComposedOp([op.adjoint() for op in reversed(self.oplist)], coeff=self.coeff)
 
-    def compose(self, other: OperatorBase) -> OperatorBase:
+    def compose(self, other: OperatorBase,
+                permutation: Optional[List[int]] = None, front: bool = False) -> OperatorBase:
+
+        new_self, other = self._expand_shorter_operator_and_permute(other, permutation)
+        new_self = cast(ComposedOp, new_self)
+
+        if front:
+            return other.compose(new_self)
         # Try composing with last element in list
         if isinstance(other, ComposedOp):
-            return ComposedOp(self.oplist + other.oplist, coeff=self.coeff * other.coeff)
+            return ComposedOp(new_self.oplist + other.oplist, coeff=new_self.coeff * other.coeff)
 
         # Try composing with last element of oplist. We only try
         # this if that last element isn't itself an
         # ComposedOp, so we can tell whether composing the
         # two elements directly worked. If it doesn't,
         # continue to the final return statement below, appending other to the oplist.
-        if not isinstance(self.oplist[-1], ComposedOp):
-            comp_with_last = self.oplist[-1].compose(other)
+        if not isinstance(new_self.oplist[-1], ComposedOp):
+            comp_with_last = new_self.oplist[-1].compose(other)
             # Attempt successful
             if not isinstance(comp_with_last, ComposedOp):
-                new_oplist = self.oplist[0:-1] + [comp_with_last]
-                return ComposedOp(new_oplist, coeff=self.coeff)
+                new_oplist = new_self.oplist[0:-1] + [comp_with_last]
+                return ComposedOp(new_oplist, coeff=new_self.coeff)
 
-        return ComposedOp(self.oplist + [other], coeff=self.coeff)
+        return ComposedOp(new_self.oplist + [other], coeff=new_self.coeff)
 
     def eval(self,
              front: Union[str, dict, np.ndarray,
diff --git a/qiskit/aqua/operators/list_ops/list_op.py b/qiskit/aqua/operators/list_ops/list_op.py
index 39ff4d278117..370d887f3e0c 100644
--- a/qiskit/aqua/operators/list_ops/list_op.py
+++ b/qiskit/aqua/operators/list_ops/list_op.py
@@ -13,15 +13,18 @@
 """ ListOp Operator Class """
 
 from functools import reduce
-from typing import List, Union, Optional, Callable, Iterator, Set, Dict
+from typing import List, Union, Optional, Callable, Iterator, Set, Dict, cast
 from numbers import Number
 
 import numpy as np
 from scipy.sparse import spmatrix
 
-from qiskit.circuit import ParameterExpression
+from qiskit.circuit import QuantumCircuit, ParameterExpression
+
 from ..legacy.base_operator import LegacyBaseOperator
 from ..operator_base import OperatorBase
+from ... import AquaError
+from ...utils import arithmetic
 
 
 class ListOp(OperatorBase):
@@ -212,11 +215,61 @@ def tensorpower(self, other: int) -> Union[OperatorBase, int]:
         from .tensored_op import TensoredOp
         return TensoredOp([self] * other)
 
-    def compose(self, other: OperatorBase) -> OperatorBase:
+    def _expand_dim(self, num_qubits: int) -> 'ListOp':
+        return ListOp([op._expand_dim(num_qubits + self.num_qubits - op.num_qubits)
+                       for op in self.oplist], combo_fn=self.combo_fn, coeff=self.coeff)
+
+    def permute(self, permutation: List[int]) -> 'ListOp':
+        """Permute the qubits of the operator.
+
+        Args:
+            permutation: A list defining where each qubit should be permuted. The qubit at index
+                j should be permuted to position permutation[j].
+
+        Returns:
+            A new ListOp representing the permuted operator.
+
+        Raises:
+            AquaError: if indices do not define a new index for each qubit.
+        """
+        new_self = self
+        circuit_size = max(permutation) + 1
+
+        try:
+            if self.num_qubits != len(permutation):
+                raise AquaError("New index must be defined for each qubit of the operator.")
+        except ValueError:
+            raise AquaError("Permute is only possible if all operators in the ListOp have the "
+                            "same number of qubits.") from ValueError
+        if self.num_qubits < circuit_size:
+            # pad the operator with identities
+            new_self = self._expand_dim(circuit_size - self.num_qubits)
+        qc = QuantumCircuit(circuit_size)
+        # extend the indices to match the size of the circuit
+        permutation \
+            = list(filter(lambda x: x not in permutation, range(circuit_size))) + permutation
+
+        # decompose permutation into sequence of transpositions
+        transpositions = arithmetic.transpositions(permutation)
+        for trans in transpositions:
+            qc.swap(trans[0], trans[1])
+
+        from qiskit.aqua.operators import CircuitOp
+
+        return CircuitOp(qc.reverse_ops()) @ new_self @ CircuitOp(qc)
+
+    def compose(self, other: OperatorBase,
+                permutation: Optional[List[int]] = None, front: bool = False) -> OperatorBase:
+
+        new_self, other = self._expand_shorter_operator_and_permute(other, permutation)
+        new_self = cast(ListOp, new_self)
+
+        if front:
+            return other.compose(new_self)
         # Avoid circular dependency
         # pylint: disable=cyclic-import,import-outside-toplevel
         from .composed_op import ComposedOp
-        return ComposedOp([self, other])
+        return ComposedOp([new_self, other])
 
     def power(self, exponent: int) -> OperatorBase:
         if not isinstance(exponent, int) or exponent <= 0:
@@ -258,8 +311,8 @@ def to_spmatrix(self) -> Union[spmatrix, List[spmatrix]]:
             [op.to_spmatrix() for op in self.oplist]) * self.coeff  # type: ignore
 
     def eval(self,
-             front: Optional[Union[str, Dict[str, complex], 'OperatorBase']] = None
-             ) -> Union['OperatorBase', float, complex, list]:
+             front: Optional[Union[str, Dict[str, complex], OperatorBase]] = None
+             ) -> Union[OperatorBase, float, complex, list]:
         """
         Evaluate the Operator's underlying function, either on a binary string or another Operator.
         A square binary Operator can be defined as a function taking a binary function to another
diff --git a/qiskit/aqua/operators/list_ops/tensored_op.py b/qiskit/aqua/operators/list_ops/tensored_op.py
index b392c0efb759..5aaab68953c6 100644
--- a/qiskit/aqua/operators/list_ops/tensored_op.py
+++ b/qiskit/aqua/operators/list_ops/tensored_op.py
@@ -54,6 +54,16 @@ def num_qubits(self) -> int:
     def distributive(self) -> bool:
         return False
 
+    def _expand_dim(self, num_qubits: int) -> 'TensoredOp':
+        """Appends I ^ num_qubits to ``oplist``. Choice of PauliOp as
+        identity is arbitrary and can be substituted for other PrimitiveOp identity.
+
+        Returns:
+            TensoredOp expanded with identity operator.
+        """
+        from qiskit.aqua.operators import I
+        return TensoredOp(self.oplist + [I ^ num_qubits], coeff=self.coeff)
+
     def tensor(self, other: OperatorBase) -> OperatorBase:
         if isinstance(other, TensoredOp):
             return TensoredOp(self.oplist + other.oplist, coeff=self.coeff * other.coeff)
diff --git a/qiskit/aqua/operators/operator_base.py b/qiskit/aqua/operators/operator_base.py
index 005d161fbe02..2571c15a92bc 100644
--- a/qiskit/aqua/operators/operator_base.py
+++ b/qiskit/aqua/operators/operator_base.py
@@ -12,7 +12,7 @@
 
 """ OperatorBase Class """
 
-from typing import Set, Union, Dict, Optional, List, cast
+from typing import Set, Union, Dict, Optional, List, cast, Tuple
 from numbers import Number
 from abc import ABC, abstractmethod
 import numpy as np
@@ -452,6 +452,32 @@ def assign_parameters(self,
         """
         raise NotImplementedError
 
+    @abstractmethod
+    def _expand_dim(self, num_qubits: int) -> 'OperatorBase':
+        """Expands the operator with identity operator of dimension 2**num_qubits.
+
+        Returns:
+            Operator corresponding to self.tensor(identity_operator), where dimension of identity
+            operator is 2 ** num_qubits.
+        """
+        raise NotImplementedError
+
+    @abstractmethod
+    def permute(self, permutation: List[int]) -> 'OperatorBase':
+        """Permutes the qubits of the operator.
+
+        Args:
+            permutation: A list defining where each qubit should be permuted. The qubit at index
+                j should be permuted to position permutation[j].
+
+        Returns:
+            A new OperatorBase containing the permuted operator.
+
+        Raises:
+            AquaError: if indices do not define a new index for each qubit.
+        """
+        raise NotImplementedError
+
     def bind_parameters(self,
                         param_dict: Dict[ParameterExpression,
                                          Union[Number,
@@ -501,6 +527,24 @@ def _get_param_dict_for_index(unrolled_dict: Dict[ParameterExpression, List[Numb
         """ Gets a single non-list-nested param_dict for a given list index from a nested one. """
         return {k: v[i] for (k, v) in unrolled_dict.items()}
 
+    def _expand_shorter_operator_and_permute(self, other: 'OperatorBase',
+                                             permutation: Optional[List[int]] = None) \
+            -> Tuple['OperatorBase', 'OperatorBase']:
+        if permutation is not None:
+            other = other.permute(permutation)
+        new_self = self
+        if not self.num_qubits == other.num_qubits:
+            # pylint: disable=cyclic-import,import-outside-toplevel
+            from .operator_globals import Zero
+            if other == Zero:
+                # Zero is special - we'll expand it to the correct qubit number.
+                other = Zero.__class__('0' * self.num_qubits)
+            elif other.num_qubits < self.num_qubits:
+                other = other._expand_dim(self.num_qubits - other.num_qubits)
+            elif other.num_qubits > self.num_qubits:
+                new_self = self._expand_dim(other.num_qubits - self.num_qubits)
+        return new_self, other
+
     # Composition
 
     def __matmul__(self, other: 'OperatorBase') -> 'OperatorBase':
@@ -515,7 +559,8 @@ def __matmul__(self, other: 'OperatorBase') -> 'OperatorBase':
         return self.compose(other)
 
     @abstractmethod
-    def compose(self, other: 'OperatorBase') -> 'OperatorBase':
+    def compose(self, other: 'OperatorBase',
+                permutation: Optional[List[int]] = None, front: bool = False) -> 'OperatorBase':
         r""" Return Operator Composition between self and other (linear algebra-style:
         A@B(x) = A(B(x))), overloaded by ``@``.
 
@@ -529,6 +574,8 @@ def compose(self, other: 'OperatorBase') -> 'OperatorBase':
 
         Args:
             other: The ``OperatorBase`` with which to compose self.
+            permutation: ``List[int]`` which defines permutation on other operator.
+            front: If front==True, return ``other.compose(self)``.
 
         Returns:
             An ``OperatorBase`` equivalent to the function composition of self and other.
diff --git a/qiskit/aqua/operators/primitive_ops/circuit_op.py b/qiskit/aqua/operators/primitive_ops/circuit_op.py
index 7d5140ccce2a..5e600d6ce472 100644
--- a/qiskit/aqua/operators/primitive_ops/circuit_op.py
+++ b/qiskit/aqua/operators/primitive_ops/circuit_op.py
@@ -108,30 +108,37 @@ def tensor(self, other: OperatorBase) -> OperatorBase:
 
         return TensoredOp([self, other])
 
-    def compose(self, other: OperatorBase) -> OperatorBase:
-        other = self._check_zero_for_composition_and_expand(other)
+    def compose(self, other: OperatorBase,
+                permutation: Optional[List[int]] = None, front: bool = False) -> OperatorBase:
+
+        new_self, other = self._expand_shorter_operator_and_permute(other, permutation)
+        new_self = cast(CircuitOp, new_self)
+
+        if front:
+            return other.compose(new_self)
+        # ignore
         # pylint: disable=cyclic-import,import-outside-toplevel
         from ..operator_globals import Zero
         from ..state_fns import CircuitStateFn
         from .pauli_op import PauliOp
         from .matrix_op import MatrixOp
 
-        if other == Zero ^ self.num_qubits:
-            return CircuitStateFn(self.primitive, coeff=self.coeff)
+        if other == Zero ^ new_self.num_qubits:
+            return CircuitStateFn(new_self.primitive, coeff=new_self.coeff)
 
         if isinstance(other, (PauliOp, CircuitOp, MatrixOp)):
             other = other.to_circuit_op()
 
         if isinstance(other, (CircuitOp, CircuitStateFn)):
-            new_qc = other.primitive.compose(self.primitive)  # type: ignore
+            new_qc = other.primitive.compose(new_self.primitive)  # type: ignore
             if isinstance(other, CircuitStateFn):
                 return CircuitStateFn(new_qc,
                                       is_measurement=other.is_measurement,
-                                      coeff=self.coeff * other.coeff)
+                                      coeff=new_self.coeff * other.coeff)
             else:
-                return CircuitOp(new_qc, coeff=self.coeff * other.coeff)
+                return CircuitOp(new_qc, coeff=new_self.coeff * other.coeff)
 
-        return super().compose(other)
+        return super(CircuitOp, new_self).compose(other)
 
     def to_matrix(self, massive: bool = False) -> np.ndarray:
         if self.num_qubits > 16 and not massive:
@@ -219,6 +226,9 @@ def reduce(self) -> OperatorBase:
                     del self.primitive.data[i]  # type: ignore
         return self
 
+    def _expand_dim(self, num_qubits: int) -> 'CircuitOp':
+        return self.permute(list(range(num_qubits, num_qubits + self.num_qubits)))
+
     def permute(self, permutation: List[int]) -> 'CircuitOp':
         r"""
         Permute the qubits of the circuit.
@@ -230,5 +240,5 @@ def permute(self, permutation: List[int]) -> 'CircuitOp':
         Returns:
             A new CircuitOp containing the permuted circuit.
         """
-        new_qc = QuantumCircuit(self.num_qubits).compose(self.primitive, qubits=permutation)
+        new_qc = QuantumCircuit(max(permutation) + 1).compose(self.primitive, qubits=permutation)
         return CircuitOp(new_qc, coeff=self.coeff)
diff --git a/qiskit/aqua/operators/primitive_ops/matrix_op.py b/qiskit/aqua/operators/primitive_ops/matrix_op.py
index 720296d9c36e..ec91a450cdce 100644
--- a/qiskit/aqua/operators/primitive_ops/matrix_op.py
+++ b/qiskit/aqua/operators/primitive_ops/matrix_op.py
@@ -12,11 +12,12 @@
 
 """ MatrixOp Class """
 
-from typing import Union, Optional, Set, Dict
+from typing import Union, Optional, Set, Dict, List, cast
 import logging
 import numpy as np
 from scipy.sparse import spmatrix
 
+from qiskit import QuantumCircuit
 from qiskit.quantum_info import Operator
 from qiskit.circuit import ParameterExpression, Instruction
 from qiskit.extensions.hamiltonian_gate import HamiltonianGate
@@ -27,6 +28,8 @@
 from ..list_ops.tensored_op import TensoredOp
 from .primitive_op import PrimitiveOp
 from ..legacy.matrix_operator import MatrixOperator
+from ...utils import arithmetic
+from ... import AquaError
 
 logger = logging.getLogger(__name__)
 
@@ -105,6 +108,10 @@ def equals(self, other: OperatorBase) -> bool:
             return self.coeff == other.coeff and self.primitive == other.primitive
         return self.coeff * self.primitive == other.coeff * other.primitive  # type: ignore
 
+    def _expand_dim(self, num_qubits: int) -> 'MatrixOp':
+        identity = np.identity(2**num_qubits, dtype=complex)
+        return MatrixOp(self.primitive.tensor(Operator(identity)), coeff=self.coeff)  # type: ignore
+
     def tensor(self, other: OperatorBase) -> OperatorBase:
         if isinstance(other.primitive, Operator):  # type: ignore
             return MatrixOp(self.primitive.tensor(other.primitive),  # type: ignore
@@ -112,14 +119,53 @@ def tensor(self, other: OperatorBase) -> OperatorBase:
 
         return TensoredOp([self, other])
 
-    def compose(self, other: OperatorBase) -> OperatorBase:
-        other = self._check_zero_for_composition_and_expand(other)
+    def compose(self, other: OperatorBase,
+                permutation: Optional[List[int]] = None, front: bool = False) -> OperatorBase:
+
+        new_self, other = self._expand_shorter_operator_and_permute(other, permutation)
+        new_self = cast(MatrixOp, new_self)
 
+        if front:
+            return other.compose(new_self)
         if isinstance(other, MatrixOp):
-            return MatrixOp(self.primitive.compose(other.primitive, front=True),  # type: ignore
-                            coeff=self.coeff * other.coeff)
+            return MatrixOp(new_self.primitive.compose(other.primitive, front=True),  # type: ignore
+                            coeff=new_self.coeff * other.coeff)
+
+        return super(MatrixOp, new_self).compose(other)
 
-        return super().compose(other)
+    def permute(self, permutation: Optional[List[int]] = None) -> 'MatrixOp':
+        """Creates a new MatrixOp that acts on the permuted qubits.
+
+        Args:
+            permutation: A list defining where each qubit should be permuted. The qubit at index
+                j should be permuted to position permutation[j].
+
+        Returns:
+            A new MatrixOp representing the permuted operator.
+
+        Raises:
+            AquaError: if indices do not define a new index for each qubit.
+        """
+        new_self = self
+        new_matrix_size = max(permutation) + 1
+
+        if self.num_qubits != len(permutation):
+            raise AquaError("New index must be defined for each qubit of the operator.")
+        if self.num_qubits < new_matrix_size:
+            # pad the operator with identities
+            new_self = self._expand_dim(new_matrix_size - self.num_qubits)
+        qc = QuantumCircuit(new_matrix_size)
+
+        # extend the indices to match the size of the new matrix
+        permutation \
+            = list(filter(lambda x: x not in permutation, range(new_matrix_size))) + permutation
+
+        # decompose permutation into sequence of transpositions
+        transpositions = arithmetic.transpositions(permutation)
+        for trans in transpositions:
+            qc.swap(trans[0], trans[1])
+        matrix = CircuitOp(qc).to_matrix()
+        return MatrixOp(matrix.transpose()) @ new_self @ MatrixOp(matrix)  # type: ignore
 
     def to_matrix(self, massive: bool = False) -> np.ndarray:
         return self.primitive.data * self.coeff  # type: ignore
diff --git a/qiskit/aqua/operators/primitive_ops/pauli_op.py b/qiskit/aqua/operators/primitive_ops/pauli_op.py
index 38edf71b2ba8..7e4adf7af6ee 100644
--- a/qiskit/aqua/operators/primitive_ops/pauli_op.py
+++ b/qiskit/aqua/operators/primitive_ops/pauli_op.py
@@ -12,7 +12,7 @@
 
 """ PauliOp Class """
 
-from typing import Union, Set, Dict, Optional, cast
+from typing import Union, Set, Dict, cast, List, Optional
 import logging
 import numpy as np
 from scipy.sparse import spmatrix
@@ -27,6 +27,7 @@
 from ..list_ops.summed_op import SummedOp
 from ..list_ops.tensored_op import TensoredOp
 from ..legacy.weighted_pauli_operator import WeightedPauliOperator
+from ... import AquaError
 
 logger = logging.getLogger(__name__)
 PAULI_GATE_MAPPING = {'X': XGate(), 'Y': YGate(), 'Z': ZGate(), 'I': IGate()}
@@ -80,6 +81,9 @@ def equals(self, other: OperatorBase) -> bool:
 
         return self.primitive == other.primitive
 
+    def _expand_dim(self, num_qubits: int) -> 'PauliOp':
+        return PauliOp(Pauli(label='I'*num_qubits).kron(self.primitive), coeff=self.coeff)
+
     def tensor(self, other: OperatorBase) -> OperatorBase:
         # Both Paulis
         if isinstance(other, PauliOp):
@@ -95,25 +99,53 @@ def tensor(self, other: OperatorBase) -> OperatorBase:
 
         return TensoredOp([self, other])
 
-    def compose(self, other: OperatorBase) -> OperatorBase:
-        other = self._check_zero_for_composition_and_expand(other)
+    def permute(self, permutation: List[int]) -> 'PauliOp':
+        """Permutes the sequence of Pauli matrices.
+
+        Args:
+            permutation: A list defining where each Pauli should be permuted. The Pauli at index
+                j of the primitive should be permuted to position permutation[j].
 
+        Returns:
+              A new PauliOp representing the permuted operator. For operator (X ^ Y ^ Z) and
+              indices=[1,2,4], it returns (X ^ I ^ Y ^ Z ^ I).
+
+        Raises:
+            AquaError: if indices do not define a new index for each qubit.
+        """
+        pauli_string = self.primitive.__str__()
+        length = max(permutation) + 1  # size of list must be +1 larger then its max index
+        new_pauli_list = ['I'] * length
+        if len(permutation) != self.num_qubits:
+            raise AquaError("List of indices to permute must have the same size as Pauli Operator")
+        for i, index in enumerate(permutation):
+            new_pauli_list[-index - 1] = pauli_string[-i - 1]
+        return PauliOp(Pauli(label=''.join(new_pauli_list)), self.coeff)
+
+    def compose(self, other: OperatorBase,
+                permutation: Optional[List[int]] = None, front: bool = False) -> OperatorBase:
+
+        new_self, other = self._expand_shorter_operator_and_permute(other, permutation)
+        new_self = cast(PauliOp, new_self)
+
+        if front:
+            return other.compose(new_self)
         # If self is identity, just return other.
-        if not any(self.primitive.x + self.primitive.z):  # type: ignore
-            return other * self.coeff  # type: ignore
+        if not any(new_self.primitive.x + new_self.primitive.z):  # type: ignore
+            return other * new_self.coeff  # type: ignore
 
         # Both Paulis
         if isinstance(other, PauliOp):
-            product, phase = Pauli.sgn_prod(self.primitive, other.primitive)
-            return PrimitiveOp(product, coeff=self.coeff * other.coeff * phase)
+            product, phase = Pauli.sgn_prod(new_self.primitive, other.primitive)
+            return PrimitiveOp(product, coeff=new_self.coeff * other.coeff * phase)
 
         # pylint: disable=cyclic-import,import-outside-toplevel
         from .circuit_op import CircuitOp
         from ..state_fns.circuit_state_fn import CircuitStateFn
         if isinstance(other, (CircuitOp, CircuitStateFn)):
-            return self.to_circuit_op().compose(other)
+            return new_self.to_circuit_op().compose(other)
 
-        return super().compose(other)
+        return super(PauliOp, new_self).compose(other)
 
     def to_matrix(self, massive: bool = False) -> np.ndarray:
         if self.num_qubits > 16 and not massive:
diff --git a/qiskit/aqua/operators/primitive_ops/primitive_op.py b/qiskit/aqua/operators/primitive_ops/primitive_op.py
index c3c1bb47af07..62adbb7a9150 100644
--- a/qiskit/aqua/operators/primitive_ops/primitive_op.py
+++ b/qiskit/aqua/operators/primitive_ops/primitive_op.py
@@ -153,29 +153,19 @@ def tensorpower(self, other: int) -> Union[OperatorBase, int]:
             temp = temp.tensor(self)
         return temp
 
-    def compose(self, other: OperatorBase) -> OperatorBase:
+    def compose(self, other: OperatorBase,
+                permutation: Optional[List[int]] = None, front: bool = False) -> \
+            OperatorBase:
         from ..list_ops.composed_op import ComposedOp
+        new_self, other = self._expand_shorter_operator_and_permute(other, permutation)
         if isinstance(other, ComposedOp):
-            comp_with_first = self.compose(other.oplist[0])
+            comp_with_first = new_self.compose(other.oplist[0])
             if not isinstance(comp_with_first, ComposedOp):
                 new_oplist = [comp_with_first] + other.oplist[1:]
                 return ComposedOp(new_oplist, coeff=other.coeff)
-            return ComposedOp([self] + other.oplist, coeff=other.coeff)  # type: ignore
-
-        return ComposedOp([self, other])
-
-    def _check_zero_for_composition_and_expand(self, other: OperatorBase) -> OperatorBase:
-        if not self.num_qubits == other.num_qubits:
-            # pylint: disable=cyclic-import,import-outside-toplevel
-            from ..operator_globals import Zero
-            if other == Zero:
-                # Zero is special - we'll expand it to the correct qubit number.
-                other = Zero.__class__('0' * self.num_qubits)
-            else:
-                raise ValueError(
-                    'Composition is not defined over Operators of different dimensions, {} and {}, '
-                    'respectively.'.format(self.num_qubits, other.num_qubits))
-        return other
+            return ComposedOp([new_self] + other.oplist, coeff=other.coeff)  # type: ignore
+
+        return ComposedOp([new_self, other])
 
     def power(self, exponent: int) -> OperatorBase:
         if not isinstance(exponent, int) or exponent <= 0:
@@ -185,6 +175,12 @@ def power(self, exponent: int) -> OperatorBase:
             temp = temp.compose(self)
         return temp
 
+    def _expand_dim(self, num_qubits: int) -> OperatorBase:
+        raise NotImplementedError
+
+    def permute(self, permutation: List[int]) -> OperatorBase:
+        raise NotImplementedError
+
     def exp_i(self) -> OperatorBase:
         """ Return Operator exponentiation, equaling e^(-i * op)"""
         # pylint: disable=cyclic-import,import-outside-toplevel
diff --git a/qiskit/aqua/operators/state_fns/circuit_state_fn.py b/qiskit/aqua/operators/state_fns/circuit_state_fn.py
index 7894bbf1b1c1..b19632e25586 100644
--- a/qiskit/aqua/operators/state_fns/circuit_state_fn.py
+++ b/qiskit/aqua/operators/state_fns/circuit_state_fn.py
@@ -137,12 +137,16 @@ def adjoint(self) -> OperatorBase:
                               coeff=np.conj(self.coeff),
                               is_measurement=(not self.is_measurement))
 
-    def compose(self, other: OperatorBase) -> OperatorBase:
-        if not self.is_measurement:
+    def compose(self, other: OperatorBase,
+                permutation: Optional[List[int]] = None, front: bool = False) -> OperatorBase:
+        if not self.is_measurement and not front:
             raise ValueError(
                 'Composition with a Statefunctions in the first operand is not defined.')
+        # type: ignore
+        new_self, other = self._expand_shorter_operator_and_permute(other, permutation)
 
-        new_self, other = self._check_zero_for_composition_and_expand(other)
+        if front:
+            return other.compose(new_self)
 
         # pylint: disable=cyclic-import,import-outside-toplevel
         from ..primitive_ops.circuit_op import CircuitOp
@@ -191,10 +195,12 @@ def tensor(self, other: OperatorBase) -> OperatorBase:
         if isinstance(other, CircuitStateFn) and other.is_measurement == self.is_measurement:
             # Avoid reimplementing tensor, just use CircuitOp's
             from ..primitive_ops.circuit_op import CircuitOp
-            from ..operator_globals import Zero
             c_op_self = CircuitOp(self.primitive, self.coeff)
             c_op_other = CircuitOp(other.primitive, other.coeff)
-            return c_op_self.tensor(c_op_other).compose(Zero)
+            c_op = c_op_self.tensor(c_op_other)
+            if isinstance(c_op, CircuitOp):
+                return CircuitStateFn(primitive=c_op.primitive, coeff=c_op.coeff,
+                                      is_measurement=self.is_measurement)
         # pylint: disable=cyclic-import
         from ..list_ops.tensored_op import TensoredOp
         return TensoredOp([self, other])
@@ -357,6 +363,11 @@ def reduce(self) -> OperatorBase:
                     del self.primitive.data[i]
         return self
 
+    def _expand_dim(self, num_qubits: int) -> 'CircuitStateFn':
+        # this is equivalent to self.tensor(identity_operator), but optimized for better performance
+        # just like in tensor method, qiskit endianness is reversed here
+        return self.permute(list(range(num_qubits, num_qubits + self.num_qubits)))
+
     def permute(self, permutation: List[int]) -> 'CircuitStateFn':
         r"""
         Permute the qubits of the circuit.
@@ -368,5 +379,5 @@ def permute(self, permutation: List[int]) -> 'CircuitStateFn':
         Returns:
             A new CircuitStateFn containing the permuted circuit.
         """
-        new_qc = QuantumCircuit(self.num_qubits).compose(self.primitive, qubits=permutation)
+        new_qc = QuantumCircuit(max(permutation) + 1).compose(self.primitive, qubits=permutation)
         return CircuitStateFn(new_qc, coeff=self.coeff, is_measurement=self.is_measurement)
diff --git a/qiskit/aqua/operators/state_fns/dict_state_fn.py b/qiskit/aqua/operators/state_fns/dict_state_fn.py
index 9e529e778a08..88e1fca8f112 100644
--- a/qiskit/aqua/operators/state_fns/dict_state_fn.py
+++ b/qiskit/aqua/operators/state_fns/dict_state_fn.py
@@ -12,18 +12,19 @@
 
 """ DictStateFn Class """
 
-from typing import Optional, Union, Set, Dict, cast
+from typing import Optional, Union, Set, Dict, cast, List
 import itertools
 import numpy as np
 from scipy import sparse
 
 from qiskit.result import Result
 from qiskit.circuit import ParameterExpression
-from qiskit.aqua import aqua_globals
+from qiskit.aqua import aqua_globals, AquaError
 
 from ..operator_base import OperatorBase
 from .state_fn import StateFn
 from ..list_ops.list_op import ListOp
+from .vector_state_fn import VectorStateFn
 
 
 class DictStateFn(StateFn):
@@ -107,6 +108,26 @@ def adjoint(self) -> OperatorBase:
                            coeff=np.conj(self.coeff),
                            is_measurement=(not self.is_measurement))
 
+    def permute(self, permutation: List[int]) -> 'DictStateFn':
+        new_num_qubits = max(permutation) + 1
+        if self.num_qubits != len(permutation):
+            raise AquaError("New index must be defined for each qubit of the operator.")
+
+        # helper function to permute the key
+        def perm(key):
+            list_key = ['0'] * new_num_qubits
+            for i, k in enumerate(permutation):
+                list_key[k] = key[i]
+            return ''.join(list_key)
+
+        new_dict = {perm(key): value for key, value in self.primitive.items()}
+        return DictStateFn(new_dict, coeff=self.coeff, is_measurement=self.is_measurement)
+
+    def _expand_dim(self, num_qubits: int) -> 'DictStateFn':
+        pad = '0'*num_qubits
+        new_dict = {key + pad: value for key, value in self.primitive.items()}
+        return DictStateFn(new_dict, coeff=self.coeff, is_measurement=self.is_measurement)
+
     def tensor(self, other: OperatorBase) -> OperatorBase:
         # Both dicts
         if isinstance(other, DictStateFn):
@@ -214,7 +235,6 @@ def eval(self,
 
         # All remaining possibilities only apply when self.is_measurement is True
 
-        from .vector_state_fn import VectorStateFn
         if isinstance(front, VectorStateFn):
             # TODO does it need to be this way for measurement?
             # return sum([v * front.primitive.data[int(b, 2)] *
diff --git a/qiskit/aqua/operators/state_fns/operator_state_fn.py b/qiskit/aqua/operators/state_fns/operator_state_fn.py
index dcd1e1b2c9b2..dcfa0a12b873 100644
--- a/qiskit/aqua/operators/state_fns/operator_state_fn.py
+++ b/qiskit/aqua/operators/state_fns/operator_state_fn.py
@@ -12,7 +12,7 @@
 
 """ OperatorStateFn Class """
 
-from typing import Union, Set
+from typing import Union, Set, List
 import numpy as np
 
 from qiskit.circuit import ParameterExpression
@@ -81,6 +81,16 @@ def adjoint(self) -> OperatorBase:
                                coeff=np.conj(self.coeff),
                                is_measurement=(not self.is_measurement))
 
+    def _expand_dim(self, num_qubits: int) -> 'OperatorStateFn':
+        return OperatorStateFn(self.primitive._expand_dim(num_qubits),
+                               coeff=self.coeff,
+                               is_measurement=self.is_measurement)
+
+    def permute(self, permutation: List[int]) -> 'OperatorStateFn':
+        return OperatorStateFn(self.primitive.permute(permutation),
+                               coeff=self.coeff,
+                               is_measurement=self.is_measurement)
+
     def tensor(self, other: OperatorBase) -> OperatorBase:
         if isinstance(other, OperatorStateFn):
             return StateFn(self.primitive.tensor(other.primitive),
diff --git a/qiskit/aqua/operators/state_fns/state_fn.py b/qiskit/aqua/operators/state_fns/state_fn.py
index e9c71a9854f2..b8dc4064ff06 100644
--- a/qiskit/aqua/operators/state_fns/state_fn.py
+++ b/qiskit/aqua/operators/state_fns/state_fn.py
@@ -12,7 +12,7 @@
 
 """ StateFn Class """
 
-from typing import Union, Optional, Callable, Set, Dict, Tuple
+from typing import Union, Optional, Callable, Set, Dict, Tuple, List
 import numpy as np
 
 from qiskit.quantum_info import Statevector
@@ -140,6 +140,21 @@ def add(self, other: OperatorBase) -> OperatorBase:
     def adjoint(self) -> OperatorBase:
         raise NotImplementedError
 
+    def _expand_dim(self, num_qubits: int) -> 'StateFn':
+        raise NotImplementedError
+
+    def permute(self, permutation: List[int]) -> OperatorBase:
+        """Permute the qubits of the state function.
+
+        Args:
+            permutation: A list defining where each qubit should be permuted. The qubit at index
+                j of the circuit should be permuted to position permutation[j].
+
+        Returns:
+            A new StateFn containing the permuted primitive.
+        """
+        raise NotImplementedError
+
     def equals(self, other: OperatorBase) -> bool:
         if not isinstance(other, type(self)) or not self.coeff == other.coeff:
             return False
@@ -188,24 +203,20 @@ def tensorpower(self, other: int) -> Union[OperatorBase, int]:
             temp = temp.tensor(self)
         return temp
 
-    def _check_zero_for_composition_and_expand(self, other: OperatorBase) \
+    def _expand_shorter_operator_and_permute(self, other: OperatorBase,
+                                             permutation: Optional[List[int]] = None) \
             -> Tuple[OperatorBase, OperatorBase]:
-        new_self = self
-        # pylint: disable=import-outside-toplevel
-        if not self.num_qubits == other.num_qubits:
-            from qiskit.aqua.operators import Zero
-            if self == StateFn({'0': 1}, is_measurement=True):
-                # Zero is special - we'll expand it to the correct qubit number.
-                new_self = StateFn('0' * self.num_qubits, is_measurement=True)
-            elif other == Zero:
-                # Zero is special - we'll expand it to the correct qubit number.
-                other = StateFn('0' * self.num_qubits)
-            else:
-                raise ValueError(
-                    'Composition is not defined over Operators of different dimensions, {} and {}, '
-                    'respectively.'.format(self.num_qubits, other.num_qubits))
-
-        return new_self, other
+
+        from qiskit.aqua.operators import Zero
+
+        if self == StateFn({'0': 1}, is_measurement=True):
+            # Zero is special - we'll expand it to the correct qubit number.
+            return StateFn('0' * other.num_qubits, is_measurement=True), other
+        elif other == Zero:
+            # Zero is special - we'll expand it to the correct qubit number.
+            return self, StateFn('0' * self.num_qubits)
+
+        return super()._expand_shorter_operator_and_permute(other, permutation)
 
     def to_matrix(self, massive: bool = False) -> np.ndarray:
         raise NotImplementedError
@@ -226,13 +237,16 @@ def to_density_matrix(self, massive: bool = False) -> np.ndarray:
         """
         raise NotImplementedError
 
-    def compose(self, other: OperatorBase) -> OperatorBase:
+    def compose(self, other: OperatorBase,
+                permutation: Optional[List[int]] = None, front: bool = False) -> OperatorBase:
         r"""
         Composition (Linear algebra-style: A@B(x) = A(B(x))) is not well defined for states
         in the binary function model, but is well defined for measurements.
 
         Args:
             other: The Operator to compose with self.
+            permutation: ``List[int]`` which defines permutation on other operator.
+            front: If front==True, return ``other.compose(self)``.
 
         Returns:
             An Operator equivalent to the function composition of self and other.
@@ -241,11 +255,14 @@ def compose(self, other: OperatorBase) -> OperatorBase:
             ValueError: If self is not a measurement, it cannot be composed from the right.
         """
         # TODO maybe allow outers later to produce density operators or projectors, but not yet.
-        if not self.is_measurement:
+        if not self.is_measurement and not front:
             raise ValueError(
                 'Composition with a Statefunction in the first operand is not defined.')
 
-        new_self, other = self._check_zero_for_composition_and_expand(other)
+        new_self, other = self._expand_shorter_operator_and_permute(other, permutation)
+
+        if front:
+            return other.compose(self)
         # TODO maybe include some reduction here in the subclasses - vector and Op, op and Op, etc.
         # pylint: disable=import-outside-toplevel
         from qiskit.aqua.operators import CircuitOp
diff --git a/qiskit/aqua/operators/state_fns/vector_state_fn.py b/qiskit/aqua/operators/state_fns/vector_state_fn.py
index e0e38355f74f..134900b0b02b 100644
--- a/qiskit/aqua/operators/state_fns/vector_state_fn.py
+++ b/qiskit/aqua/operators/state_fns/vector_state_fn.py
@@ -13,8 +13,9 @@
 """ VectorStateFn Class """
 
 
-from typing import Union, Set, Optional, Dict
+from typing import Union, Set, Optional, Dict, List
 import numpy as np
+from qiskit import QuantumCircuit
 
 from qiskit.quantum_info import Statevector
 from qiskit.circuit import ParameterExpression
@@ -23,6 +24,7 @@
 from ..operator_base import OperatorBase
 from .state_fn import StateFn
 from ..list_ops.list_op import ListOp
+from ...utils import arithmetic
 
 
 class VectorStateFn(StateFn):
@@ -76,6 +78,51 @@ def adjoint(self) -> OperatorBase:
                              coeff=np.conj(self.coeff),
                              is_measurement=(not self.is_measurement))
 
+    def permute(self, permutation: List[int]) -> 'VectorStateFn':
+        new_self = self
+        new_num_qubits = max(permutation) + 1
+
+        if self.num_qubits != len(permutation):
+            # raise AquaError("New index must be defined for each qubit of the operator.")
+            pass
+        if self.num_qubits < new_num_qubits:
+            # pad the operator with identities
+            new_self = self._expand_dim(new_num_qubits - self.num_qubits)
+        qc = QuantumCircuit(new_num_qubits)
+
+        # extend the permutation indices to match the size of the new matrix
+        permutation \
+            = list(filter(lambda x: x not in permutation, range(new_num_qubits))) + permutation
+
+        # decompose permutation into sequence of transpositions
+        transpositions = arithmetic.transpositions(permutation)
+        for trans in transpositions:
+            qc.swap(trans[0], trans[1])
+        from .. import CircuitOp
+        matrix = CircuitOp(qc).to_matrix()
+        vector = new_self.primitive.data
+        return VectorStateFn(primitive=matrix.dot(vector),
+                             coeff=self.coeff,
+                             is_measurement=self.is_measurement)
+
+    def to_dict_fn(self) -> StateFn:
+        """Creates the equivalent state function of type DictStateFn.
+
+        Returns:
+            A new DictStateFn equivalent to ``self``.
+        """
+        from .dict_state_fn import DictStateFn
+
+        num_qubits = self.num_qubits
+        new_dict = {format(i, 'b').zfill(num_qubits): v for i, v in enumerate(self.primitive.data)}
+        return DictStateFn(new_dict, coeff=self.coeff, is_measurement=self.is_measurement)
+
+    def _expand_dim(self, num_qubits: int) -> 'VectorStateFn':
+        primitive = np.zeros(2**num_qubits, dtype=complex)
+        return VectorStateFn(self.primitive.tensor(primitive),
+                             coeff=self.coeff,
+                             is_measurement=self.is_measurement)
+
     def tensor(self, other: OperatorBase) -> OperatorBase:
         if isinstance(other, VectorStateFn):
             return StateFn(self.primitive.tensor(other.primitive),
@@ -145,9 +192,9 @@ def eval(self,
 
         # pylint: disable=cyclic-import,import-outside-toplevel
         from ..operator_globals import EVAL_SIG_DIGITS
-        from .dict_state_fn import DictStateFn
         from .operator_state_fn import OperatorStateFn
         from .circuit_state_fn import CircuitStateFn
+        from .dict_state_fn import DictStateFn
         if isinstance(front, DictStateFn):
             return np.round(sum([v * self.primitive.data[int(b, 2)] * front.coeff  # type: ignore
                                  for (b, v) in front.primitive.items()]) * self.coeff,
diff --git a/qiskit/aqua/utils/arithmetic.py b/qiskit/aqua/utils/arithmetic.py
index 96caf6a64f39..a4eb20d2688a 100644
--- a/qiskit/aqua/utils/arithmetic.py
+++ b/qiskit/aqua/utils/arithmetic.py
@@ -14,6 +14,7 @@
 Arithmetic Utilities
 """
 
+from typing import List, Tuple
 import numpy as np
 
 
@@ -93,3 +94,40 @@ def next_power_of_2_base(n):
         base += 1
 
     return base
+
+
+def transpositions(permutation: List[int]) -> List[Tuple[int, int]]:
+    """Return a sequence of transpositions, corresponding to the permutation.
+
+    Args:
+        permutation: The ``List[int]`` defining the permutation. An element at index ``j`` should be
+            permuted to index ``permutation[j]``.
+
+    Returns:
+        List of transpositions, corresponding to the permutation. For permutation = [3, 0, 2, 1],
+        returns [(0,1), (0,3)]
+    """
+    unchecked = [True] * len(permutation)
+    cyclic_form = []
+    for i in range(len(permutation)):
+        if unchecked[i]:
+            cycle = [i]
+            unchecked[i] = False
+            j = i
+            while unchecked[permutation[j]]:
+                j = permutation[j]
+                cycle.append(j)
+                unchecked[j] = False
+            if len(cycle) > 1:
+                cyclic_form.append(cycle)
+    cyclic_form.sort()
+    res = []
+    for x in cyclic_form:
+        len_x = len(x)
+        if len_x == 2:
+            res.append((x[0], x[1]))
+        elif len_x > 2:
+            first = x[0]
+            for y in x[len_x - 1:0:-1]:
+                res.append((first, y))
+    return res