diff --git a/src/Simulation/Native/src/simulator/capi.cpp b/src/Simulation/Native/src/simulator/capi.cpp
index 5f081acb159..344445f46c4 100644
--- a/src/Simulation/Native/src/simulator/capi.cpp
+++ b/src/Simulation/Native/src/simulator/capi.cpp
@@ -21,13 +21,13 @@ MICROSOFT_QUANTUM_DECL void destroy(_In_ unsigned id)
 
 MICROSOFT_QUANTUM_DECL void seed(_In_ unsigned id, _In_ unsigned s)
 {
-  psis[id]->seed(s);
+  Microsoft::Quantum::Simulator::get(id)->seed(s);
 }
 
 // non-quantum
 MICROSOFT_QUANTUM_DECL std::size_t random_choice(_In_ unsigned id, _In_ std::size_t n, _In_reads_(n) double* p)
 {
-    return psis[id]->random(n, p);
+    return Microsoft::Quantum::Simulator::get(id)->random(n, p);
 }
 
 MICROSOFT_QUANTUM_DECL double JointEnsembleProbability(_In_ unsigned id, _In_ unsigned n, _In_reads_(n) int* b, _In_reads_(n) unsigned* q)
@@ -36,35 +36,35 @@ MICROSOFT_QUANTUM_DECL double JointEnsembleProbability(_In_ unsigned id, _In_ un
     for (unsigned i = 0; i < n; ++i)
         bv.push_back(static_cast<Gates::Basis>(*(b + i)));
     std::vector<unsigned> qv(q, q + n);
-    return psis[id]->JointEnsembleProbability( bv, qv);
+    return Microsoft::Quantum::Simulator::get(id)->JointEnsembleProbability( bv, qv);
 }
 
 
 MICROSOFT_QUANTUM_DECL void allocateQubit(_In_ unsigned id, _In_ unsigned q)
 {
-    psis[id]->allocateQubit(q);
+    Microsoft::Quantum::Simulator::get(id)->allocateQubit(q);
 }
 
 MICROSOFT_QUANTUM_DECL void release(_In_ unsigned id, _In_ unsigned q)
 {
-    psis[id]->release(q);
+    Microsoft::Quantum::Simulator::get(id)->release(q);
 }
 
 MICROSOFT_QUANTUM_DECL unsigned num_qubits(_In_ unsigned id)
 {
-  return psis[id]->num_qubits();
+  return Microsoft::Quantum::Simulator::get(id)->num_qubits();
 }
 
 #define FWDGATE1(G)                                                                                                    \
 MICROSOFT_QUANTUM_DECL void G(_In_ unsigned id, _In_ unsigned q)                                                                          \
     {                                                                                                                  \
-        psis[id]->G(q);                                                                           \
+        Microsoft::Quantum::Simulator::get(id)->G(q);                                                                           \
     }
 #define FWDCSGATE1(G)                                                                                                  \
     MICROSOFT_QUANTUM_DECL void MC##G(_In_ unsigned id, _In_ unsigned n, _In_reads_(n) unsigned* c, _In_ unsigned q)                                             \
     {                                                                                                                  \
         std::vector<unsigned> vc(c, c + n);                                                                            \
-        psis[id]->C##G(vc, q);                                                                    \
+        Microsoft::Quantum::Simulator::get(id)->C##G(vc, q);                                                                    \
     }
 #define FWD(G) FWDGATE1(G) FWDCSGATE1(G)
 
@@ -88,14 +88,14 @@ FWD(AdjT)
     
 MICROSOFT_QUANTUM_DECL void R(_In_ unsigned id, _In_ unsigned b, _In_ double phi, _In_ unsigned q)
 {
-    psis[id]->R(static_cast<Gates::Basis>(b), phi, q);
+    Microsoft::Quantum::Simulator::get(id)->R(static_cast<Gates::Basis>(b), phi, q);
 }
   
 // multi-controlled rotations
 MICROSOFT_QUANTUM_DECL void MCR(_In_ unsigned id, _In_ unsigned b, _In_ double phi, _In_ unsigned nc, _In_reads_(nc) unsigned* c, _In_ unsigned q)
 {
     std::vector<unsigned> cv(c, c + nc);
-    psis[id]->CR(static_cast<Gates::Basis>(b), phi, cv, q);
+    Microsoft::Quantum::Simulator::get(id)->CR(static_cast<Gates::Basis>(b), phi, cv, q);
 }
 
 // Exponential of Pauli operators
@@ -105,7 +105,7 @@ MICROSOFT_QUANTUM_DECL void Exp(_In_ unsigned id, _In_ unsigned n, _In_reads_(n)
     for (unsigned i = 0; i < n; ++i)
         bv.push_back(static_cast<Gates::Basis>(*(b + i)));
     std::vector<unsigned> qv(q, q + n);
-    psis[id]->Exp(bv, phi, qv);
+    Microsoft::Quantum::Simulator::get(id)->Exp(bv, phi, qv);
 }
 MICROSOFT_QUANTUM_DECL void MCExp(_In_ unsigned id, _In_ unsigned n, _In_reads_(n) unsigned* b, _In_ double phi, _In_ unsigned nc, _In_reads_(nc) unsigned* c, _In_reads_(n) unsigned* q)
 {
@@ -114,13 +114,13 @@ MICROSOFT_QUANTUM_DECL void MCExp(_In_ unsigned id, _In_ unsigned n, _In_reads_(
         bv.push_back(static_cast<Gates::Basis>(*(b + i)));
     std::vector<unsigned> qv(q, q + n);
     std::vector<unsigned> cv(c, c + nc);
-    psis[id]->CExp(bv, phi, cv, qv);
+    Microsoft::Quantum::Simulator::get(id)->CExp(bv, phi, cv, qv);
 }
     
 // measurements
 MICROSOFT_QUANTUM_DECL unsigned M(_In_ unsigned id, _In_ unsigned q)
 {
-    return (unsigned)psis[id]->M(q);
+    return (unsigned)Microsoft::Quantum::Simulator::get(id)->M(q);
 }
 MICROSOFT_QUANTUM_DECL unsigned Measure(_In_ unsigned id, _In_ unsigned n, _In_reads_(n) unsigned* b, _In_reads_(n) unsigned* q)
 {
@@ -128,7 +128,7 @@ MICROSOFT_QUANTUM_DECL unsigned Measure(_In_ unsigned id, _In_ unsigned n, _In_r
     for (unsigned i = 0; i < n; ++i)
         bv.push_back(static_cast<Gates::Basis>(*(b + i)));
     std::vector<unsigned> qv(q, q + n);
-    return (unsigned)psis[id]->Measure(bv, qv);
+    return (unsigned)Microsoft::Quantum::Simulator::get(id)->Measure(bv, qv);
 }
     
 // apply permutation of basis states to the wave function
@@ -136,33 +136,33 @@ MICROSOFT_QUANTUM_DECL void PermuteBasis(_In_ unsigned id, _In_ unsigned n, _In_
     _In_reads_(table_size) std::size_t *permutation_table)
 {
     const std::vector<unsigned> qs(q, q + n);
-    psis[id]->permuteBasis(qs, table_size, permutation_table, false);
+    Microsoft::Quantum::Simulator::get(id)->permuteBasis(qs, table_size, permutation_table, false);
 }
 MICROSOFT_QUANTUM_DECL void AdjPermuteBasis(_In_ unsigned id, _In_ unsigned n, _In_reads_(n) unsigned* q, _In_ std::size_t table_size,
     _In_reads_(table_size) std::size_t *permutation_table)
 {
     const std::vector<unsigned> qs(q, q + n);
-    psis[id]->permuteBasis(qs, table_size, permutation_table, true);
+    Microsoft::Quantum::Simulator::get(id)->permuteBasis(qs, table_size, permutation_table, true);
 }
 
 
 // dump wavefunction to given callback until callback returns false
 MICROSOFT_QUANTUM_DECL void Dump(_In_ unsigned id, _In_ bool (*callback)(size_t, double, double))
 {
-    psis[id]->dump(callback);
+    Microsoft::Quantum::Simulator::get(id)->dump(callback);
 }
 
 // dump the wavefunction of the subset of qubits to the given callback returns false
 MICROSOFT_QUANTUM_DECL bool DumpQubits(_In_ unsigned id, _In_ unsigned n, _In_reads_(n) unsigned* q, _In_ bool(*callback)(size_t, double, double))
 {
     std::vector<unsigned> qs(q, q + n);
-    return psis[id]->dumpQubits(qs, callback);
+    return Microsoft::Quantum::Simulator::get(id)->dumpQubits(qs, callback);
 }
 
 // dump the list of logical qubit ids to given callback
 MICROSOFT_QUANTUM_DECL void DumpIds(_In_ unsigned id, _In_ void(*callback)(unsigned))
 {
-    psis[id]->dumpIds(callback);
+    Microsoft::Quantum::Simulator::get(id)->dumpIds(callback);
 }
 
 }
diff --git a/src/Simulation/Native/src/simulator/factory.cpp b/src/Simulation/Native/src/simulator/factory.cpp
index 6d46fbcf7f7..b80b02b35f6 100644
--- a/src/Simulation/Native/src/simulator/factory.cpp
+++ b/src/Simulation/Native/src/simulator/factory.cpp
@@ -5,6 +5,7 @@
 #include "config.hpp"
 #include "util/cpuid.hpp"
 #include <iostream>
+#include <shared_mutex>
 
 namespace Microsoft
 {
@@ -31,7 +32,8 @@ namespace Microsoft
   {
     namespace Simulator
     {
-      mutex_type _mutex;
+      std::shared_mutex _mutex;
+      std::vector<std::shared_ptr<SimulatorInterface>> _psis;
 
       SimulatorInterface* createSimulator(unsigned maxlocal)
       {
@@ -52,26 +54,26 @@ namespace Microsoft
 
       MICROSOFT_QUANTUM_DECL unsigned create(unsigned maxlocal)
       {
-        std::lock_guard<mutex_type> lock(_mutex);
+        std::lock_guard<std::shared_mutex> lock(_mutex);
 
         size_t emptySlot = -1;
-        for (auto const& s : psis) 
+        for (auto const& s : _psis) 
         {
             if (s == NULL) 
             {
-                emptySlot = &s - &psis[0];
+                emptySlot = &s - &_psis[0];
                 break;
             }
         }
 
         if (emptySlot == -1) 
         {
-            psis.push_back(std::shared_ptr<SimulatorInterface>(createSimulator(maxlocal)));
-            emptySlot = psis.size() - 1;
+            _psis.push_back(std::shared_ptr<SimulatorInterface>(createSimulator(maxlocal)));
+            emptySlot = _psis.size() - 1;
         }
         else 
         {
-            psis[emptySlot] = std::shared_ptr<SimulatorInterface>(createSimulator(maxlocal));
+            _psis[emptySlot] = std::shared_ptr<SimulatorInterface>(createSimulator(maxlocal));
         }
 
         return static_cast<unsigned>(emptySlot);
@@ -79,12 +81,17 @@ namespace Microsoft
 
       MICROSOFT_QUANTUM_DECL void destroy(unsigned id)
       {
-        std::lock_guard<mutex_type> lock(_mutex);
+        std::lock_guard<std::shared_mutex> lock(_mutex);
 
-        psis[id].reset();
+        _psis[id].reset();
       }
 
-      MICROSOFT_QUANTUM_DECL std::vector<std::shared_ptr<SimulatorInterface>> psis;
+      MICROSOFT_QUANTUM_DECL std::shared_ptr<SimulatorInterface>& get(unsigned id)
+      {
+        std::shared_lock<std::shared_mutex> shared_lock(_mutex);
+
+        return _psis[id];
+      }
 
     }
   }
diff --git a/src/Simulation/Native/src/simulator/factory.hpp b/src/Simulation/Native/src/simulator/factory.hpp
index ec9d097fc8a..979b4a3b052 100644
--- a/src/Simulation/Native/src/simulator/factory.hpp
+++ b/src/Simulation/Native/src/simulator/factory.hpp
@@ -12,8 +12,7 @@ namespace Microsoft
     {
       MICROSOFT_QUANTUM_DECL unsigned create(unsigned =0u);
       MICROSOFT_QUANTUM_DECL void destroy(unsigned);
-      
-      extern MICROSOFT_QUANTUM_DECL std::vector<std::shared_ptr<SimulatorInterface>> psis;
+      MICROSOFT_QUANTUM_DECL std::shared_ptr<SimulatorInterface>& get(unsigned);
     }
   }
 }
diff --git a/src/Simulation/Native/src/simulator/factory_test.cpp b/src/Simulation/Native/src/simulator/factory_test.cpp
index 6c196205713..d4745fa2d90 100644
--- a/src/Simulation/Native/src/simulator/factory_test.cpp
+++ b/src/Simulation/Native/src/simulator/factory_test.cpp
@@ -8,7 +8,7 @@ using namespace Microsoft::Quantum::Simulator;
 
 int main(int argc, char** argv)
 {
-  auto& sim = psis[create()];
+  auto sim = get(create());
 
   unsigned q=0; // qubit number
   sim->allocateQubit(q);