diff --git a/src/framework/operations.hpp b/src/framework/operations.hpp
index ea8b68da1a..484dd3660a 100644
--- a/src/framework/operations.hpp
+++ b/src/framework/operations.hpp
@@ -383,7 +383,6 @@ enum class OpType {
   superop,
   roerror,
   noise_switch,
-  sample_noise,
   // Save instructions
   save_state,
   save_expval,
@@ -532,9 +531,6 @@ inline std::ostream &operator<<(std::ostream &stream, const OpType &type) {
   case OpType::qerror_loc:
     stream << "qerror_loc";
     break;
-  case OpType::sample_noise:
-    stream << "sample_noise";
-    break;
   case OpType::noise_switch:
     stream << "noise_switch";
     break;
@@ -614,11 +610,14 @@ struct Op {
       string_params; // used for label, control-flow, and boolean functions
 
   // Conditional Operations
-  bool conditional = false; // is gate conditional gate
-  uint_t conditional_reg; // (opt) the (single) register location to look up for
-                          // conditional
-  BinaryOp binary_op;     // (opt) boolean function relation
-  std::shared_ptr<CExpr> expr; // (opt) classical expression
+  // is gate conditional gate
+  bool conditional = false;
+  // (opt) the (single) register location to look up for conditional
+  uint_t conditional_reg = 0;
+  // (opt) boolean function relation
+  BinaryOp binary_op;
+  // (opt) classical expression
+  std::shared_ptr<CExpr> expr = nullptr;
 
   // Measurement
   reg_t memory;    // (opt) register operation it acts on (measure)
@@ -640,6 +639,9 @@ struct Op {
   // Save
   DataSubType save_type = DataSubType::single;
 
+  // runtime noise sampling
+  bool sample_noise = false;
+
   // runtime parameter bind
   bool has_bind_params = false;
 };
@@ -1319,8 +1321,14 @@ inline Op bind_parameter(const Op &src, const uint_t iparam,
   op.type = src.type;
   op.name = src.name;
   op.qubits = src.qubits;
+  op.regs = src.regs;
+  op.int_params = src.int_params;
+  op.string_params = src.string_params;
   op.conditional = src.conditional;
   op.conditional_reg = src.conditional_reg;
+  op.binary_op = src.binary_op;
+  op.expr = src.expr;
+  op.has_bind_params = false;
 
   if (src.params.size() > 0) {
     uint_t stride = src.params.size() / num_params;
diff --git a/src/noise/noise_model.hpp b/src/noise/noise_model.hpp
index 23dbccc9bd..4063aaa58a 100644
--- a/src/noise/noise_model.hpp
+++ b/src/noise/noise_model.hpp
@@ -65,7 +65,13 @@ class NoiseModel {
                        const Method method = Method::circuit,
                        bool sample_at_runtime = false) const;
 
-  NoiseOps sample_noise_loc(const Operations::Op &op, RngEngine &rng) const;
+  NoiseOps sample_noise_at_runtime(const Operations::Op &op,
+                                   RngEngine &rng) const;
+
+  void sample_noise_at_runtime(const Operations::Op &op,
+                               NoiseModel::NoiseOps &noise_before,
+                               NoiseModel::NoiseOps &noise_after,
+                               RngEngine &rng) const;
 
   // Enable superop sampling method
   // This will cause all QuantumErrors stored in the noise model
@@ -151,14 +157,12 @@ class NoiseModel {
   void sample_local_quantum_noise(const Operations::Op &op,
                                   NoiseOps &noise_before, NoiseOps &noise_after,
                                   RngEngine &rng, const Method method,
-                                  const reg_t &mapping,
-                                  bool sample_at_runtime) const;
+                                  const reg_t &mapping) const;
 
   void sample_nonlocal_quantum_noise(const Operations::Op &op,
                                      NoiseOps &noise_ops, NoiseOps &noise_after,
                                      RngEngine &rng, const Method method,
-                                     const reg_t &mapping,
-                                     bool sample_at_runtime) const;
+                                     const reg_t &mapping) const;
 
   // Sample noise for the current operation
   NoiseOps sample_noise_helper(const Operations::Op &op, RngEngine &rng,
@@ -176,9 +180,6 @@ class NoiseModel {
                                   const std::vector<reg_t> &op_qubits,
                                   const std::vector<reg_t> &noise_qubits);
 
-  // create loc noise
-  NoiseOps create_noise_loc(const Operations::Op &op) const;
-
   // Flags which say whether the local or nonlocal error tables are used
   bool local_quantum_errors_ = false;
   bool nonlocal_quantum_errors_ = false;
@@ -245,7 +246,8 @@ class NoiseModel {
   std::unordered_set<Method> enabled_methods_ =
       std::unordered_set<Method>({Method::circuit});
 
-  // saved qubit mapping for runtime loc
+  // saved qubit mapping and method for runtime noise sampling
+  mutable Method circ_method_;
   mutable reg_t circ_mapping_;
 };
 
@@ -281,21 +283,57 @@ Circuit NoiseModel::sample_noise(const Circuit &circ, RngEngine &rng,
   return sample_noise_circuit(circ, rng, method, sample_at_runtime);
 }
 
-NoiseModel::NoiseOps NoiseModel::sample_noise_loc(const Operations::Op &op,
-                                                  RngEngine &rng) const {
-  auto noise_ops =
-      sample_noise_op(op, rng, Method::circuit, circ_mapping_, false);
+NoiseModel::NoiseOps
+NoiseModel::sample_noise_at_runtime(const Operations::Op &op,
+                                    RngEngine &rng) const {
+  // Return operator set
+  NoiseOps noise_before;
+  NoiseOps noise_after;
+
+  sample_noise_at_runtime(op, noise_before, noise_after, rng);
+
+  // Combine errors
+  auto &noise_ops = noise_before;
+  noise_ops.reserve(noise_before.size() + noise_after.size() + 1);
+  Operations::Op op_sampled = op;
+  op_sampled.sample_noise = false;
+  noise_ops.push_back(op_sampled);
+  noise_ops.insert(noise_ops.end(),
+                   std::make_move_iterator(noise_after.begin()),
+                   std::make_move_iterator(noise_after.end()));
+
+  return noise_ops;
+}
+
+void NoiseModel::sample_noise_at_runtime(const Operations::Op &op,
+                                         NoiseModel::NoiseOps &noise_before,
+                                         NoiseModel::NoiseOps &noise_after,
+                                         RngEngine &rng) const {
+  // Apply local errors first
+  sample_local_quantum_noise(op, noise_before, noise_after, rng, circ_method_,
+                             circ_mapping_);
+  // Apply nonlocal errors second
+  sample_nonlocal_quantum_noise(op, noise_before, noise_after, rng,
+                                circ_method_, circ_mapping_);
+  // Apply readout error to measure ops
+  if (op.type == Operations::OpType::measure) {
+    sample_readout_noise(op, noise_after, rng, circ_mapping_);
+  }
 
   // If original op is conditional, make all the noise operations also
   // conditional
   if (op.conditional) {
-    for (auto &noise_op : noise_ops) {
+    for (auto &noise_op : noise_before) {
+      noise_op.conditional = op.conditional;
+      noise_op.conditional_reg = op.conditional_reg;
+      noise_op.binary_op = op.binary_op;
+    }
+    for (auto &noise_op : noise_after) {
       noise_op.conditional = op.conditional;
       noise_op.conditional_reg = op.conditional_reg;
       noise_op.binary_op = op.binary_op;
     }
   }
-  return noise_ops;
 }
 
 Circuit NoiseModel::sample_noise_circuit(const Circuit &circ, RngEngine &rng,
@@ -355,6 +393,7 @@ Circuit NoiseModel::sample_noise_circuit(const Circuit &circ, RngEngine &rng,
   noisy_circ.set_params();
   if (sample_at_runtime) {
     noisy_circ.can_sample = false;
+    circ_method_ = method;
     circ_mapping_ = mapping;
   }
   return noisy_circ;
@@ -518,21 +557,27 @@ NoiseModel::sample_noise_helper(const Operations::Op &op, RngEngine &rng,
   NoiseOps noise_after;
   // Apply local errors first
   sample_local_quantum_noise(op, noise_before, noise_after, rng, method,
-                             mapping, sample_at_runtime);
+                             mapping);
   // Apply nonlocal errors second
   sample_nonlocal_quantum_noise(op, noise_before, noise_after, rng, method,
-                                mapping, sample_at_runtime);
+                                mapping);
   // Apply readout error to measure ops
   if (op.type == Operations::OpType::measure) {
     sample_readout_noise(op, noise_after, rng, mapping);
   }
 
+  if (sample_at_runtime &&
+      (noise_before.size() > 0 || noise_after.size() > 0)) {
+    NoiseOps ret(1);
+    ret[0] = op;
+    ret[0].sample_noise = true;
+    return ret;
+  }
+
   // Combine errors
   auto &noise_ops = noise_before;
   noise_ops.reserve(noise_before.size() + noise_after.size() + 1);
-  if (op.type != Operations::OpType::sample_noise) {
-    noise_ops.push_back(op);
-  }
+  noise_ops.push_back(op);
   noise_ops.insert(noise_ops.end(),
                    std::make_move_iterator(noise_after.begin()),
                    std::make_move_iterator(noise_after.end()));
@@ -655,8 +700,7 @@ void NoiseModel::sample_local_quantum_noise(const Operations::Op &op,
                                             NoiseOps &noise_before,
                                             NoiseOps &noise_after,
                                             RngEngine &rng, const Method method,
-                                            const reg_t &mapping,
-                                            bool sample_at_runtime) const {
+                                            const reg_t &mapping) const {
   // If no errors are defined pass
   if (local_quantum_errors_ == false)
     return;
@@ -714,11 +758,7 @@ void NoiseModel::sample_local_quantum_noise(const Operations::Op &op,
                                     : iter_default->second;
         for (auto &pos : error_positions) {
           NoiseOps noise_ops;
-          if (sample_at_runtime)
-            noise_ops = create_noise_loc(op);
-          else
-            noise_ops =
-                quantum_errors_[pos].sample_noise(op_qubits, rng, method);
+          noise_ops = quantum_errors_[pos].sample_noise(op_qubits, rng, method);
           // Duplicate same sampled error operations
           if (quantum_errors_[pos].errors_after())
             noise_after.insert(noise_after.end(), noise_ops.begin(),
@@ -734,8 +774,7 @@ void NoiseModel::sample_local_quantum_noise(const Operations::Op &op,
 
 void NoiseModel::sample_nonlocal_quantum_noise(
     const Operations::Op &op, NoiseOps &noise_before, NoiseOps &noise_after,
-    RngEngine &rng, const Method method, const reg_t &mapping,
-    bool sample_at_runtime) const {
+    RngEngine &rng, const Method method, const reg_t &mapping) const {
   // If no errors are defined pass
   if (nonlocal_quantum_errors_ == false)
     return;
@@ -783,12 +822,8 @@ void NoiseModel::sample_nonlocal_quantum_noise(
           auto &error_positions = target_pair.second;
           for (auto &pos : error_positions) {
             NoiseOps ops;
-            if (sample_at_runtime)
-              ops = create_noise_loc(op);
-            else {
-              ops = quantum_errors_[pos].sample_noise(string2reg(target_qubits),
-                                                      rng, method);
-            }
+            ops = quantum_errors_[pos].sample_noise(string2reg(target_qubits),
+                                                    rng, method);
             if (quantum_errors_[pos].errors_after())
               noise_after.insert(noise_after.end(), ops.begin(), ops.end());
             else
@@ -800,14 +835,6 @@ void NoiseModel::sample_nonlocal_quantum_noise(
   }
 }
 
-NoiseModel::NoiseOps
-NoiseModel::create_noise_loc(const Operations::Op &op) const {
-  NoiseOps ops(1);
-  ops[0] = op;
-  ops[0].type = Operations::OpType::sample_noise;
-  return ops;
-}
-
 cmatrix_t NoiseModel::op2superop(const Operations::Op &op) const {
   switch (op.type) {
   case Operations::OpType::superop:
diff --git a/src/simulators/batch_shots_executor.hpp b/src/simulators/batch_shots_executor.hpp
index 612e5ed289..8e829f753f 100644
--- a/src/simulators/batch_shots_executor.hpp
+++ b/src/simulators/batch_shots_executor.hpp
@@ -99,6 +99,9 @@ class BatchShotsExecutor : public virtual MultiStateExecutor<state_t> {
                                InputIterator last_meas, uint_t shots,
                                uint_t i_group, ResultItr result,
                                std::vector<RngEngine> &rng);
+
+  // check if ops contains pauli ops
+  bool check_pauli_only(std::vector<Operations::Op> &ops);
 };
 
 template <class state_t>
@@ -506,13 +509,13 @@ void BatchShotsExecutor<state_t>::apply_ops_batched_shots_for_group(
 #endif
 
   for (auto op = first; op != last; ++op) {
-    if (op->type == Operations::OpType::sample_noise) {
+    if (op->sample_noise) {
       if (op->expr) {
         for (uint_t j = Base::top_state_of_group_[i_group];
              j < Base::top_state_of_group_[i_group + 1]; j++) {
           Base::states_[j].qreg().enable_batch(false);
           Base::states_[j].qreg().read_measured_data(Base::states_[j].creg());
-          std::vector<Operations::Op> nops = noise.sample_noise_loc(
+          std::vector<Operations::Op> nops = noise.sample_noise_at_runtime(
               *op, rng[j - Base::top_state_of_group_[i_group]]);
           for (uint_t k = 0; k < nops.size(); k++) {
             Base::states_[j].apply_op(
@@ -526,55 +529,87 @@ void BatchShotsExecutor<state_t>::apply_ops_batched_shots_for_group(
 
       // sample error here
       uint_t count = Base::num_states_in_group_[i_group];
-      std::vector<std::vector<Operations::Op>> noise_ops(count);
+      std::vector<std::vector<Operations::Op>> noise_ops_before(count);
+      std::vector<std::vector<Operations::Op>> noise_ops_after(count);
 
-      uint_t count_ops = 0;
-      uint_t non_pauli_gate_count = 0;
+      uint_t count_ops_before = 0;
+      uint_t count_ops_after = 0;
+      bool pauli_only_before = true;
+      bool pauli_only_after = true;
       if (num_inner_threads > 1) {
-#pragma omp parallel for reduction(+: count_ops,non_pauli_gate_count) num_threads(num_inner_threads)
+#pragma omp parallel for reduction(+: count_ops_before, count_ops_after) reduction(&: pauli_only_before, pauli_only_after) num_threads(num_inner_threads)
         for (int_t j = 0; j < (int_t)count; j++) {
-          noise_ops[j] = noise.sample_noise_loc(*op, rng[j]);
-
-          if (!(noise_ops[j].size() == 0 ||
-                (noise_ops[j].size() == 1 && noise_ops[j][0].name == "id"))) {
-            count_ops++;
-            for (uint_t k = 0; k < noise_ops[j].size(); k++) {
-              if (noise_ops[j][k].name != "id" && noise_ops[j][k].name != "x" &&
-                  noise_ops[j][k].name != "y" && noise_ops[j][k].name != "z" &&
-                  noise_ops[j][k].name != "pauli") {
-                non_pauli_gate_count++;
-                break;
-              }
-            }
+          noise.sample_noise_at_runtime(*op, noise_ops_before[j],
+                                        noise_ops_after[j], rng[j]);
+
+          pauli_only_before &= check_pauli_only(noise_ops_before[j]);
+          pauli_only_after &= check_pauli_only(noise_ops_after[j]);
+
+          if (!(noise_ops_before[j].size() == 0 ||
+                (noise_ops_before[j].size() == 1 &&
+                 noise_ops_before[j][0].name == "id"))) {
+            count_ops_before++;
+          }
+          if (!(noise_ops_after[j].size() == 0 ||
+                (noise_ops_after[j].size() == 1 &&
+                 noise_ops_after[j][0].name == "id"))) {
+            count_ops_after++;
           }
         }
       } else {
         for (uint_t j = 0; j < count; j++) {
-          noise_ops[j] = noise.sample_noise_loc(*op, rng[j]);
-
-          if (!(noise_ops[j].size() == 0 ||
-                (noise_ops[j].size() == 1 && noise_ops[j][0].name == "id"))) {
-            count_ops++;
-            for (uint_t k = 0; k < noise_ops[j].size(); k++) {
-              if (noise_ops[j][k].name != "id" && noise_ops[j][k].name != "x" &&
-                  noise_ops[j][k].name != "y" && noise_ops[j][k].name != "z" &&
-                  noise_ops[j][k].name != "pauli") {
-                non_pauli_gate_count++;
-                break;
-              }
-            }
+          noise.sample_noise_at_runtime(*op, noise_ops_before[j],
+                                        noise_ops_after[j], rng[j]);
+
+          pauli_only_before &= check_pauli_only(noise_ops_before[j]);
+          pauli_only_after &= check_pauli_only(noise_ops_after[j]);
+
+          if (!(noise_ops_before[j].size() == 0 ||
+                (noise_ops_before[j].size() == 1 &&
+                 noise_ops_before[j][0].name == "id"))) {
+            count_ops_before++;
+          }
+          if (!(noise_ops_after[j].size() == 0 ||
+                (noise_ops_after[j].size() == 1 &&
+                 noise_ops_after[j][0].name == "id"))) {
+            count_ops_after++;
           }
         }
       }
 
-      if (count_ops == 0) {
-        continue; // do nothing
+      // noise before op
+      if (count_ops_before > 0) {
+        if (pauli_only_before) { // optimization for Pauli error
+          Base::states_[istate].qreg().apply_batched_pauli_ops(
+              noise_ops_before);
+        } else {
+          // otherwise execute each circuit
+          apply_batched_noise_ops(i_group, noise_ops_before, result_it, rng);
+        }
       }
-      if (non_pauli_gate_count == 0) { // optimization for Pauli error
-        Base::states_[istate].qreg().apply_batched_pauli_ops(noise_ops);
-      } else {
-        // otherwise execute each circuit
-        apply_batched_noise_ops(i_group, noise_ops, result_it, rng);
+      // apply original op
+      if (op->expr || !apply_batched_op(istate, *op, result_it, rng,
+                                        final_ops && (op + 1 == last))) {
+        // call apply_op for each state
+        for (uint_t j = 0; j < Base::num_states_in_group_[i_group]; j++) {
+          uint_t is = Base::top_state_of_group_[i_group] + j;
+          uint_t ip = (Base::global_state_index_ + is) /
+                      Base::num_shots_per_bind_param_;
+          Base::states_[is].qreg().enable_batch(false);
+          Base::states_[is].qreg().read_measured_data(Base::states_[is].creg());
+          Base::states_[is].apply_op(*op, *(result_it + ip), rng[j],
+                                     final_ops && (op + 1 == last));
+          Base::states_[is].qreg().enable_batch(true);
+        }
+      }
+      // noise after op
+      if (count_ops_after > 0) {
+        if (pauli_only_after) { // optimization for Pauli error
+          Base::states_[istate].qreg().apply_batched_pauli_ops(noise_ops_after);
+        } else {
+          // otherwise execute each circuit
+          apply_batched_noise_ops(i_group, noise_ops_after, result_it, rng);
+        }
       }
     } else {
       if (!op->expr && apply_batched_op(istate, *op, result_it, rng,
@@ -596,6 +631,18 @@ void BatchShotsExecutor<state_t>::apply_ops_batched_shots_for_group(
   }
 }
 
+template <class state_t>
+bool BatchShotsExecutor<state_t>::check_pauli_only(
+    std::vector<Operations::Op> &ops) {
+  for (uint_t k = 0; k < ops.size(); k++) {
+    if (ops[k].name != "id" && ops[k].name != "x" && ops[k].name != "y" &&
+        ops[k].name != "z" && ops[k].name != "pauli") {
+      return false;
+    }
+  }
+  return true;
+}
+
 template <class state_t>
 void BatchShotsExecutor<state_t>::apply_batched_noise_ops(
     const int_t i_group, const std::vector<std::vector<Operations::Op>> &ops,
diff --git a/src/simulators/multi_state_executor.hpp b/src/simulators/multi_state_executor.hpp
index 3f180bea1b..7f80ae2645 100644
--- a/src/simulators/multi_state_executor.hpp
+++ b/src/simulators/multi_state_executor.hpp
@@ -468,7 +468,7 @@ void MultiStateExecutor<state_t>::run_circuit_with_shot_branching(
   // initial state
   waiting_branches.push_back(std::make_shared<Branch>());
   waiting_branches[0]->set_shots(shots_storage);
-  waiting_branches[0]->op_iterator() = first;
+  waiting_branches[0]->set_iterator(first);
   if (Base::num_bind_params_ > 1) {
     waiting_branches[0]->set_param_index(global_state_index_ + ishot,
                                          Base::num_shots_per_bind_param_);
@@ -494,7 +494,7 @@ void MultiStateExecutor<state_t>::run_circuit_with_shot_branching(
         break;
       uint_t sid = top_state + i;
       waiting_branches[i]->state_index() = sid;
-      waiting_branches[i]->op_iterator() = first;
+      waiting_branches[i]->set_iterator(first);
       branches.push_back(waiting_branches[i]);
 
       // initialize state
@@ -514,7 +514,7 @@ void MultiStateExecutor<state_t>::run_circuit_with_shot_branching(
     while (num_active_states > 0) { // loop until all branches execute all ops
       // functor for ops execution
       auto apply_ops_func = [this, &branches, &noise, &par_results, measure_seq,
-                             last, par_shots, num_active_states](int_t i) {
+                             last, par_shots](int_t i) {
         uint_t istate, state_end;
         istate = branches.size() * i / par_shots;
         state_end = branches.size() * (i + 1) / par_shots;
@@ -524,56 +524,9 @@ void MultiStateExecutor<state_t>::run_circuit_with_shot_branching(
         for (; istate < state_end; istate++) {
           state_t &state = states_[branches[istate]->state_index()];
 
-          while (branches[istate]->op_iterator() != measure_seq ||
-                 branches[istate]->additional_ops().size() > 0) {
-            // execute additional ops first if avaiable
-            if (branches[istate]->additional_ops().size() > 0) {
-              int_t iadd = 0;
-              int_t num_add = branches[istate]->additional_ops().size();
-              while (iadd < num_add) {
-                if (apply_branching_op(*branches[istate],
-                                       branches[istate]->additional_ops()[iadd],
-                                       par_results[i].begin(), false)) {
-                  // check if there are new branches
-                  if (branches[istate]->num_branches() > 0) {
-                    // if there are additional ops remaining, queue them on new
-                    // branches
-                    for (uint_t k = iadd + 1;
-                         k < branches[istate]->additional_ops().size(); k++) {
-                      for (uint_t l = 0; l < branches[istate]->num_branches();
-                           l++)
-                        branches[istate]->branches()[l]->add_op_after_branch(
-                            branches[istate]->additional_ops()[k]);
-                    }
-                    branches[istate]->remove_empty_branches();
-                    state.creg() = branches[istate]->creg();
-
-                    // if there are some branches still remaining
-                    if (branches[istate]->num_branches() > 0) {
-                      nbranch += branches[istate]->num_branches();
-                      break;
-                    }
-                    iadd = 0;
-                    num_add = branches[istate]->additional_ops().size();
-                  }
-                } else {
-                  state.apply_op(branches[istate]->additional_ops()[iadd],
-                                 par_results[i][0], dummy_rng, false);
-                }
-                iadd++;
-              }
-              branches[istate]->clear_additional_ops();
-              // if there are some branches still remaining
-              if (branches[istate]->num_branches() > 0) {
-                nbranch += branches[istate]->num_branches();
-                break;
-              }
-            }
+          while (branches[istate]->op_iterator() != measure_seq) {
             OpItr op = branches[istate]->op_iterator();
-            if (op == measure_seq)
-              break;
 
-            // then execute ops
             if (!state.creg().check_conditional(*op)) {
               branches[istate]->advance_iterator();
               continue;
@@ -581,13 +534,13 @@ void MultiStateExecutor<state_t>::run_circuit_with_shot_branching(
             if (branches[istate]->apply_control_flow(state.creg(), measure_seq))
               continue;
 
-            // runtime noise sampling
-            if (op->type == Operations::OpType::sample_noise) {
+            branches[istate]->advance_iterator();
+            if (op->sample_noise) {
+              // runtime noise sampling
               branches[istate]->apply_runtime_noise_sampling(state.creg(), *op,
                                                              noise);
-            }
-            // runtime parameterizaion
-            else if (op->has_bind_params) {
+            } else if (op->has_bind_params) {
+              // runtime parameterizaion
               apply_runtime_parameterization(*branches[istate], *op);
             } else {
               if (!apply_branching_op(*branches[istate], *op,
@@ -598,7 +551,6 @@ void MultiStateExecutor<state_t>::run_circuit_with_shot_branching(
               }
             }
 
-            branches[istate]->advance_iterator();
             if (branches[istate]->num_branches() > 0) {
               branches[istate]->remove_empty_branches();
               state.creg() = branches[istate]->creg();
@@ -672,8 +624,7 @@ void MultiStateExecutor<state_t>::run_circuit_with_shot_branching(
       // check if there are remaining ops
       num_active_states = 0;
       for (uint_t i = 0; i < branches.size(); i++) {
-        if (branches[i]->op_iterator() != measure_seq ||
-            branches[i]->additional_ops().size() > 0)
+        if (branches[i]->op_iterator() != measure_seq)
           num_active_states++;
       }
     }
diff --git a/src/simulators/shot_branching.hpp b/src/simulators/shot_branching.hpp
index 0d81f707a4..ab6805ec25 100644
--- a/src/simulators/shot_branching.hpp
+++ b/src/simulators/shot_branching.hpp
@@ -39,6 +39,7 @@ class Branch {
 
   // additional operations applied after shot branching
   std::vector<Operations::Op> additional_ops_;
+  uint_t additional_op_pos_;
 
   // mark for control flow
   std::unordered_map<std::string, OpItr> flow_marks_;
@@ -49,8 +50,11 @@ class Branch {
   // branches from this
   std::vector<std::shared_ptr<Branch>> branches_;
 
+  // this flag is used for initialize op
+  bool initialize_after_reset_ = false;
+
 public:
-  Branch(void) {}
+  Branch(void) { additional_op_pos_ = 0; }
   ~Branch() {
     shots_.clear();
     additional_ops_.clear();
@@ -61,22 +65,28 @@ class Branch {
     creg_ = src.creg_;
     iter_ = src.iter_;
     flow_marks_ = src.flow_marks_;
+    additional_ops_ = src.additional_ops_;
+    additional_op_pos_ = src.additional_op_pos_;
   }
 
   uint_t &state_index(void) { return state_index_; }
   uint_t &root_state_index(void) { return root_state_index_; }
   ClassicalRegister &creg(void) { return creg_; }
   std::vector<RngEngine> &rng_shots(void) { return shots_; }
-  OpItr &op_iterator(void) { return iter_; }
   std::unordered_map<std::string, OpItr> &marks(void) { return flow_marks_; }
   uint_t num_branches(void) { return branches_.size(); }
   std::vector<std::shared_ptr<Branch>> &branches(void) { return branches_; }
+  bool &initialize_after_reset(void) { return initialize_after_reset_; }
+
+  void set_iterator(OpItr &iter) { iter_ = iter; }
+  OpItr op_iterator(void);
 
   uint_t num_shots(void) { return shots_.size(); }
   void clear(void) {
     shots_.clear();
     additional_ops_.clear();
     branches_.clear();
+    additional_op_pos_ = 0;
   }
   void clear_branch(void) { branches_.clear(); }
 
@@ -91,16 +101,16 @@ class Branch {
   void add_op_after_branch(Operations::Op &op) {
     additional_ops_.push_back(op);
   }
-  void copy_ops_after_branch(std::vector<Operations::Op> &ops) {
-    additional_ops_ = ops;
+  void add_ops_after_branch(std::vector<Operations::Op> &ops) {
+    additional_ops_.insert(additional_ops_.end(), ops.begin(), ops.end());
   }
-  void clear_additional_ops(void) { additional_ops_.clear(); }
-
-  std::vector<Operations::Op> &additional_ops(void) { return additional_ops_; }
 
   void branch_shots(reg_t &shots, int_t nbranch);
 
   bool apply_control_flow(ClassicalRegister &creg, OpItr last) {
+    if (additional_ops_.size() > additional_op_pos_)
+      return false;
+
     if (iter_->type == Operations::OpType::mark) {
       flow_marks_[iter_->string_params[0]] = iter_;
       iter_++;
@@ -236,26 +246,35 @@ void Branch::branch_shots_by_params(void) {
   }
 }
 
-void Branch::advance_iterator(void) {
-  iter_++;
-  for (uint_t i = 0; i < branches_.size(); i++) {
-    branches_[i]->iter_++;
+OpItr Branch::op_iterator(void) {
+  if (additional_ops_.size() > additional_op_pos_) {
+    OpItr it = additional_ops_.cbegin();
+    it += additional_op_pos_;
+    return it;
   }
+  return iter_;
+}
+
+void Branch::advance_iterator(void) {
+  if (additional_ops_.size() > additional_op_pos_)
+    additional_op_pos_++;
+  else
+    iter_++;
 }
 
 bool Branch::apply_runtime_noise_sampling(const ClassicalRegister &creg,
                                           const Operations::Op &op,
                                           const Noise::NoiseModel &noise) {
-  if (op.type != Operations::OpType::sample_noise)
-    return false;
-
   uint_t nshots = num_shots();
   reg_t shot_map(nshots);
   std::vector<std::vector<Operations::Op>> noises;
 
+  if (!op.sample_noise)
+    return false;
+
   for (uint_t i = 0; i < nshots; i++) {
     std::vector<Operations::Op> noise_ops =
-        noise.sample_noise_loc(op, shots_[i]);
+        noise.sample_noise_at_runtime(op, shots_[i]);
 
     // search same noise ops
     int_t pos = -1;
@@ -264,6 +283,11 @@ bool Branch::apply_runtime_noise_sampling(const ClassicalRegister &creg,
         continue;
       bool same = true;
       for (uint_t k = 0; k < noise_ops.size(); k++) {
+        if (noise_ops[k].sample_noise) {
+          noise_ops[k].sample_noise = false;
+          continue; // skip original op
+        }
+
         if (noise_ops[k].type != noises[j][k].type ||
             noise_ops[k].name != noises[j][k].name)
           same = false;
@@ -323,6 +347,14 @@ bool Branch::apply_runtime_noise_sampling(const ClassicalRegister &creg,
       }
     }
 
+    if (noises.size() == 0) {
+      for (uint_t k = 0; k < noise_ops.size(); k++) {
+        if (noise_ops[k].sample_noise) {
+          noise_ops[k].sample_noise = false;
+        }
+      }
+    }
+
     if (pos < 0) { // if not found, add noise ops to the list
       shot_map[i] = noises.size();
       noises.push_back(noise_ops);
@@ -334,37 +366,51 @@ bool Branch::apply_runtime_noise_sampling(const ClassicalRegister &creg,
   creg_ = creg;
   branch_shots(shot_map, noises.size());
   for (uint_t i = 0; i < noises.size(); i++) {
-    branches_[i]->copy_ops_after_branch(noises[i]);
+    branches_[i]->add_ops_after_branch(noises[i]);
   }
 
   return true;
 }
 
 void Branch::remove_empty_branches(void) {
-  int_t istart = 0;
+  // find first branch that has at least one shot
+  int_t iroot = -1;
   for (uint_t j = 0; j < branches_.size(); j++) {
     if (branches_[j]->num_shots() > 0) {
-      // copy shots to the root
-      shots_ = branches_[j]->rng_shots();
-      param_index_ = branches_[j]->param_index_;
-      param_shots_ = branches_[j]->param_shots_;
-      additional_ops_ = branches_[j]->additional_ops();
-      creg_ = branches_[j]->creg();
-      branches_[j].reset();
-      istart = j + 1;
+      iroot = j;
       break;
     }
     branches_[j].reset();
   }
 
-  std::vector<std::shared_ptr<Branch>> new_branches;
+  // copy shots to the root
+  shots_ = branches_[iroot]->rng_shots();
+  param_index_ = branches_[iroot]->param_index_;
+  param_shots_ = branches_[iroot]->param_shots_;
+  creg_ = branches_[iroot]->creg();
+  initialize_after_reset_ = branches_[iroot]->initialize_after_reset_;
 
-  for (uint_t j = istart; j < branches_.size(); j++) {
-    if (branches_[j]->num_shots() > 0)
-      new_branches.push_back(branches_[j]);
-    else
+  std::vector<std::shared_ptr<Branch>> new_branches;
+  for (uint_t j = iroot; j < branches_.size(); j++) {
+    if (branches_[j]->num_shots() > 0) {
+      // update additional ops if there are remaining additional ops
+      if (additional_ops_.size() > additional_op_pos_) {
+        branches_[j]->additional_ops_.insert(
+            branches_[j]->additional_ops_.end(),
+            additional_ops_.begin() + additional_op_pos_,
+            additional_ops_.end());
+      }
+      if (j != iroot)
+        new_branches.push_back(branches_[j]);
+    } else
       branches_[j].reset();
   }
+
+  additional_ops_ = branches_[iroot]->additional_ops_;
+  additional_op_pos_ = 0;
+  branches_[iroot].reset();
+  branches_.clear();
+
   branches_ = new_branches;
 }
 
diff --git a/src/simulators/statevector/statevector_executor.hpp b/src/simulators/statevector/statevector_executor.hpp
index 826642fc6a..26b00e801e 100644
--- a/src/simulators/statevector/statevector_executor.hpp
+++ b/src/simulators/statevector/statevector_executor.hpp
@@ -1647,19 +1647,22 @@ void Executor<state_t>::apply_initialize(CircuitExecutor::Branch &root,
     }
   }
 
-  if (root.additional_ops().size() == 0) {
+  if (!root.initialize_after_reset()) {
     apply_reset(root, qubits);
 
-    Operations::Op op;
-    op.type = OpType::initialize;
-    op.name = "initialize";
-    op.qubits = qubits;
-    op.params = params;
-    for (uint_t i = 0; i < root.num_branches(); i++) {
-      root.branches()[i]->add_op_after_branch(op);
+    if (root.num_branches() > 0) {
+      Operations::Op op;
+      op.type = OpType::initialize;
+      op.name = "initialize";
+      op.qubits = qubits;
+      op.params = params;
+      for (uint_t i = 0; i < root.num_branches(); i++) {
+        root.branches()[i]->add_op_after_branch(op);
+        root.branches()[i]->initialize_after_reset() = true;
+      }
+      return; // initialization will be done in next call because of shot
+              // branching in reset
     }
-    return; // initialization will be done in next call because of shot
-            // branching in reset
   }
 
   Base::states_[root.state_index()].qreg().initialize_component(qubits, params);
diff --git a/src/simulators/tensor_network/tensor_net_executor.hpp b/src/simulators/tensor_network/tensor_net_executor.hpp
index 5bcc47532f..cd147e0797 100644
--- a/src/simulators/tensor_network/tensor_net_executor.hpp
+++ b/src/simulators/tensor_network/tensor_net_executor.hpp
@@ -276,19 +276,22 @@ void Executor<state_t>::apply_initialize(CircuitExecutor::Branch &root,
     }
   }
 
-  if (root.additional_ops().size() == 0) {
+  if (!root.initialize_after_reset()) {
     apply_reset(root, qubits);
 
-    Operations::Op op;
-    op.type = OpType::initialize;
-    op.name = "initialize";
-    op.qubits = qubits;
-    op.params = params;
-    for (uint_t i = 0; i < root.num_branches(); i++) {
-      root.branches()[i]->add_op_after_branch(op);
+    if (root.num_branches() > 0) {
+      Operations::Op op;
+      op.type = OpType::initialize;
+      op.name = "initialize";
+      op.qubits = qubits;
+      op.params = params;
+      for (uint_t i = 0; i < root.num_branches(); i++) {
+        root.branches()[i]->add_op_after_branch(op);
+        root.branches()[i]->initialize_after_reset() = true;
+      }
+      return; // initialization will be done in next call because of shot
+              // branching in reset
     }
-    return; // initialization will be done in next call because of shot
-            // branching in reset
   }
 
   Base::states_[root.state_index()].qreg().initialize_component(qubits, params);
diff --git a/src/transpile/fusion.hpp b/src/transpile/fusion.hpp
index c65160fed9..9948c640dc 100644
--- a/src/transpile/fusion.hpp
+++ b/src/transpile/fusion.hpp
@@ -174,7 +174,7 @@ class UnitaryFusion : public FusionMethod {
   };
 
   virtual bool can_apply(const op_t &op, uint_t max_fused_qubits) const {
-    if (op.conditional)
+    if (op.conditional || op.sample_noise)
       return false;
     switch (op.type) {
     case optype_t::matrix:
@@ -221,7 +221,7 @@ class SuperOpFusion : public UnitaryFusion {
   };
 
   virtual bool can_apply(const op_t &op, uint_t max_fused_qubits) const {
-    if (op.conditional)
+    if (op.conditional || op.sample_noise)
       return false;
     switch (op.type) {
     case optype_t::kraus:
@@ -271,7 +271,7 @@ class KrausFusion : public UnitaryFusion {
   };
 
   virtual bool can_apply(const op_t &op, uint_t max_fused_qubits) const {
-    if (op.conditional)
+    if (op.conditional || op.sample_noise)
       return false;
     switch (op.type) {
     case optype_t::kraus: