diff --git a/qiskit/aqua/operators/primitive_ops/circuit_op.py b/qiskit/aqua/operators/primitive_ops/circuit_op.py
index da7006d22f..7d5140ccce 100644
--- a/qiskit/aqua/operators/primitive_ops/circuit_op.py
+++ b/qiskit/aqua/operators/primitive_ops/circuit_op.py
@@ -16,7 +16,8 @@
 import logging
 import numpy as np
 
-from qiskit import QuantumCircuit, BasicAer, execute
+import qiskit
+from qiskit import QuantumCircuit
 from qiskit.circuit.library import IGate
 from qiskit.circuit import Instruction, ParameterExpression
 
@@ -139,14 +140,7 @@ def to_matrix(self, massive: bool = False) -> np.ndarray:
                 ' in this case {0}x{0} elements.'
                 ' Set massive=True if you want to proceed.'.format(2 ** self.num_qubits))
 
-        # NOTE: not reversing qubits!! We generally reverse endianness when converting between
-        # circuit or Pauli representation and matrix representation, but we don't need to here
-        # because the Unitary simulator already presents the endianness of the circuit unitary in
-        # forward endianness.
-        unitary_backend = BasicAer.get_backend('unitary_simulator')
-        unitary = execute(self.to_circuit(),
-                          unitary_backend,
-                          optimization_level=0).result().get_unitary()
+        unitary = qiskit.quantum_info.Operator(self.to_circuit()).data
         # pylint: disable=cyclic-import
         from ..operator_globals import EVAL_SIG_DIGITS
         return np.round(unitary * self.coeff, decimals=EVAL_SIG_DIGITS)
diff --git a/qiskit/aqua/operators/primitive_ops/pauli_op.py b/qiskit/aqua/operators/primitive_ops/pauli_op.py
index 4910a9b851..38edf71b2b 100644
--- a/qiskit/aqua/operators/primitive_ops/pauli_op.py
+++ b/qiskit/aqua/operators/primitive_ops/pauli_op.py
@@ -220,13 +220,13 @@ def exp_i(self) -> OperatorBase:
                 else self.coeff
             # Y rotation
             if corrected_x[sig_qubit_index] and corrected_z[sig_qubit_index]:
-                rot_op = PrimitiveOp(RYGate(coeff))
+                rot_op = PrimitiveOp(RYGate(2 * coeff))
             # Z rotation
             elif corrected_z[sig_qubit_index]:
-                rot_op = PrimitiveOp(RZGate(coeff))
+                rot_op = PrimitiveOp(RZGate(2 * coeff))
             # X rotation
             elif corrected_x[sig_qubit_index]:
-                rot_op = PrimitiveOp(RXGate(coeff))
+                rot_op = PrimitiveOp(RXGate(2 * coeff))
 
             from ..operator_globals import I
             left_pad = I.tensorpower(sig_qubit_index)
diff --git a/releasenotes/notes/pauli-expi-c0c97395190622c0.yaml b/releasenotes/notes/pauli-expi-c0c97395190622c0.yaml
new file mode 100644
index 0000000000..7eaed2bd8d
--- /dev/null
+++ b/releasenotes/notes/pauli-expi-c0c97395190622c0.yaml
@@ -0,0 +1,9 @@
+---
+fixes:
+  - |
+    Make PauliOp.exp_i() generate the correct matrix with the following changes.
+    1) There was previously an error in the phase of a factor of 2.
+    2) The global phase was ignored when converting the circuit
+    to a matrix. We now use qiskit.quantum_info.Operator, which is
+    generally useful for converting a circuit to a unitary matrix,
+    when possible.
diff --git a/test/aqua/operators/test_evolution.py b/test/aqua/operators/test_evolution.py
index 85d019ffba..fbd7e52b1c 100644
--- a/test/aqua/operators/test_evolution.py
+++ b/test/aqua/operators/test_evolution.py
@@ -18,6 +18,7 @@
 import numpy as np
 import scipy.linalg
 
+import qiskit
 from qiskit.circuit import ParameterVector, Parameter
 
 from qiskit.aqua.operators import (X, Y, Z, I, CX, H, ListOp, CircuitOp, Zero, EvolutionFactory,
@@ -29,6 +30,13 @@
 class TestEvolution(QiskitAquaTestCase):
     """Evolution tests."""
 
+    def test_exp_i(self):
+        """ exponential of Pauli test """
+        op = Z.exp_i()
+        gate = op.to_circuit().data[0][0]
+        self.assertIsInstance(gate, qiskit.circuit.library.RZGate)
+        self.assertEqual(gate.params[0], 2)
+
     def test_pauli_evolution(self):
         """ pauli evolution test """
         op = (-1.052373245772859 * I ^ I) + \
diff --git a/test/aqua/operators/test_op_construction.py b/test/aqua/operators/test_op_construction.py
index a3850cef18..d12de2ba92 100644
--- a/test/aqua/operators/test_op_construction.py
+++ b/test/aqua/operators/test_op_construction.py
@@ -17,6 +17,7 @@
 
 from test.aqua import QiskitAquaTestCase
 import itertools
+import scipy
 from scipy.stats import unitary_group
 import numpy as np
 from ddt import ddt, data
@@ -201,6 +202,14 @@ def test_to_matrix(self):
         np.testing.assert_array_almost_equal(
             op6.to_matrix(), op5.to_matrix() + Operator.from_label('+r').data)
 
+    def test_circuit_op_to_matrix(self):
+        """ test CircuitOp.to_matrix """
+        qc = QuantumCircuit(1)
+        qc.rz(1.0, 0)
+        qcop = CircuitOp(qc)
+        np.testing.assert_array_almost_equal(
+            qcop.to_matrix(), scipy.linalg.expm(-0.5j * Z.to_matrix()))
+
     def test_matrix_to_instruction(self):
         """Test MatrixOp.to_instruction yields an Instruction object."""
         matop = (H ^ 3).to_matrix_op()